"""
核反应堆子通道分析程序 - 基于两流体模型

本程序实现了一个简化的核反应堆子通道分析工具，用于模拟和分析燃料棒束中冷却剂的流动特性。
采用两流体模型处理气液两相流动，考虑了轴向和横向的质量、动量传递。

主要功能：
1. 模拟子通道内的两相流动
2. 计算空泡率、压力分布
3. 分析液相和气相的速度场
4. 考虑子通道间的横向混合

物理模型：
- 两流体模型：分别处理气相和液相
- 考虑相间作用力
- 包含重力和压力梯度效应
- 简化的横向混合模型

数值方法：
- 采用显式时间推进格式
- 交错网格离散
- 迭代求解直至收敛
- 使用稀疏矩阵求解线性方程组
作者：[李泉伯 2024-12-12]
"""

import numpy as np
from CoolProp.CoolProp import PropsSI
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import os
import csv
from matplotlib.animation import FuncAnimation
import multiprocessing as mp
from functools import partial
from scipy.sparse import csr_matrix, lil_matrix
from scipy.sparse.linalg import gmres
import numpy.linalg as la
from numba import jit, njit, prange
import numba as nb
import json
import traceback
import scipy.optimize as opt
import logging
import sys


# 设置matplotlib的中文显示
plt.rcParams['font.sans-serif'] = ['SimHei']  
plt.rcParams['axes.unicode_minus'] = False  


class SubchannelAnalysis:
    """
    子通道分析的主类，实现了两流体模型的求解
    
    该类实现了基于两流体模型的核反应堆子通道分析，可以模拟和分析3×3棒束16子通道
    系统中的冷却剂流动特性。
    
    主要属性：
        n_axial (int): 轴向节点数
        n_channels (int): 子通道数量
        channel_connections (dict): 子通道连接关系
        connection_coefficients (dict): 子通道间的连接系数
        
    主要方法：
        solve_mass_momentum(): 求解质量和动量守恒方程
        calculate_cross_flow(): 计算横向流动
        get_fluid_properties(): 获取流体物性
        run_simulation(): 运行完整的模拟计算
    """
    
    def __init__(self, n_axial=20):
        """
        初始化3×3棒束16子通道分析器
        
        参数：
            n_axial (int): 轴向网格节点数，默认为2000
            
        初始化的主要内容：
            1. 几何参数设置
            2. 工况参数定义
            3. 网格划分
            4. 场变量初始化
            5. 子通道连接关系建立
        """
        # === 基本几何参数（单位：米）===
        self.rod_diameter = 14.5e-3  # 燃料棒直径
        self.rod_pitch = 18.7e-3    # 燃料棒间距
        self.bundle_size = 3        # 3×3棒束
        self.n_channels = 16        # 总子通道数
        self.n_axial = n_axial      # 轴向网格数
        self.channel_length = 4.0   # 通道总长度 (m)
        self.dz = self.channel_length / n_axial  # 轴向网格间距 (m)
        
        # === 入口工况参数 ===
        self.inlet_pressure = 15.5e6    # 入口压力 (Pa)
        self.inlet_temp = 583.15        # 入口温度 (K)
        self.inlet_mass_flux = 3500.0   # 入口质量流量 (kg/m²·s)
        self.outlet_pressure = 15.3e6   # 出口压力 (Pa)
        self.average_heat_flux = 1.2e6  # 平均热流密度 (W/m²)
        
        # === 时间步长相关参数 ===
        self.dt = 1e-4  # 时间步长 (s)
        self.current_time = 0.0  # 当前时间 (s)
        self.save_interval = 100  # 保存间隔
        
        # === 定义子通道类型（16子通道布局）===
        self.channel_types = {
            1: 'corner',   # 左上角
            2: 'edge',     # 上边
            3: 'edge',     # 上边
            4: 'corner',   # 右上角
            5: 'edge',     # 左边
            6: 'center',   # 中心
            7: 'center',   # 中心
            8: 'edge',     # 右边
            9: 'edge',     # 左边
            10: 'center',  # 中心
            11: 'center',  # 中心
            12: 'edge',    # 右边
            13: 'corner',  # 左下角
            14: 'edge',    # 下边
            15: 'edge',    # 下边
            16: 'corner'   # 右下角
        }
        
        # === 子通道几何参数定义 ===
        self.channel_geometry = {
            'corner': {
                'flow_area': 5.050e-5,        # 流通面积 (m²)
                'wetted_perimeter': 2.835e-2,  # 湿周长 (m)
                'heated_perimeter': 1.14e-2    # 加热周长 (m)
            },
            'edge': {
                'flow_area': 1.177e-4,        # 流通面积 (m²)
                'wetted_perimeter': 4.149e-2,  # 湿周长 (m)
                'heated_perimeter': 2.28e-2    # 加热周长 (m)
            },
            'center': {
                'flow_area': 1.868e-4,        # 流通面积 (m²)
                'wetted_perimeter': 4.548e-2,  # 湿周长 (m)
                'heated_perimeter': 4.56e-2    # 加热周长 (m)
            }
        }
        
        # === 计算总流通面积和质量流量 ===
        self.total_flow_area = sum(self.channel_geometry[self.channel_types[j+1]]['flow_area'] 
                                  for j in range(self.n_channels))
        
        # === 定义子通道连接关系（16子通道拓扑结构）===
        self.channel_connections = {
            1: [2, 5],           # 1号通道连接2,5
            2: [1, 3, 6],        # 2号通道连接1,3,6
            3: [2, 4, 7],        # 3号通道连接2,4,7
            4: [3, 8],           # 4号通道连接3,8
            5: [1, 6, 9],        # 5号通道连接1,6,9
            6: [2, 5, 7, 10],    # 6号通道连接2,5,7,10
            7: [3, 6, 8, 11],    # 7号通道连接3,6,8,11
            8: [4, 7, 12],       # 8号通道连接4,7,12
            9: [5, 10, 13],      # 9号通道连接5,10,13
            10: [6, 9, 11, 14],  # 10号通道连接6,9,11,14
            11: [7, 10, 12, 15], # 11号通道连接7,10,12,15
            12: [8, 11, 16],     # 12号通道连接8,11,16
            13: [9, 14],         # 13号通道连接9,14
            14: [13, 15, 10],    # 14号通道连接13,15,10
            15: [14, 16, 11],    # 15号通道连接14,16,11
            16: [15, 12]         # 16号通道连接15,12
        }
        
        # === 计算连接系数（gap宽度）===
        self.connection_gaps = {}
        for i in range(1, 17):
            for j in self.channel_connections[i]:
                if (i, j) not in self.connection_gaps and (j, i) not in self.connection_gaps:
                    # 根据子通道类型确定连接宽度
                    if self.channel_types[i] == 'corner' or self.channel_types[j] == 'corner':
                        gap_width = 0.02  # 角通道连接宽度 (m)
                    else:
                        gap_width = 0.025 # 其他通道连接宽度 (m)
                    self.connection_gaps[(i, j)] = gap_width
                    self.connection_gaps[(j, i)] = gap_width
        
        # === 创建结果保存文件夹 ===
        self.results_folder = 'transient_results'
        if not os.path.exists(self.results_folder):
            os.makedirs(self.results_folder)
        
        # === 初始化场变量数组 ===
        self.initialize_field_variables()
        
        # === 入口工况参数设置 ===
        self.inlet_pressure = 15.5e6  # 入口压力 (Pa)
        self.inlet_temp = 300         # 入口温度 (K)
        self.mass_flux = 1500         # 质量流量 (kg/m²s)
        
        # === 场变量初始化 ===
        # 空泡率场初始化
        self.void_fraction = np.zeros((n_axial, self.n_channels))    # 空泡率 (-)
        self.void_fraction[0, :] = 0.1  # 入口处初始空泡率设为10%

        # 速度场初始化
        self.liquid_velocity = np.zeros((n_axial, self.n_channels))  # 液相速度 (m/s)
        self.vapor_velocity = np.zeros((n_axial, self.n_channels))   # 气相速度 (m/s)

        # 根据质量流量计算入口速度
        rho_l = self.get_fluid_properties(self.inlet_pressure, self.inlet_temp, 'liquid')['density']
        self.liquid_velocity[0, :] = self.mass_flux / rho_l * (1 - self.void_fraction[0, :])  
        self.vapor_velocity[0, :] = self.liquid_velocity[0, :] * 0.1  # 气相速度初始为液相的10%

        # 压力和温度场初始化
        self.pressure = np.full((n_axial, self.n_channels), self.inlet_pressure)  # 压力场 (Pa)
        self.temperature = np.full((n_axial, self.n_channels), self.inlet_temp)   # 温度场 (K)
        self.liquid_enthalpy = np.zeros((n_axial, self.n_channels))  # 液相焓 (J/kg)
        self.vapor_enthalpy = np.zeros((n_axial, self.n_channels))   # 气相焓 (J/kg)
        
        # === 两流体模型物理量初始化 ===
        # 相分数
        self.alpha_l = np.zeros((n_axial, self.n_channels))  # 液相分数 (-)
        self.alpha_v = np.zeros((n_axial, self.n_channels))  # 气相分数 (-)
        
        # 速度场
        self.u_l = np.zeros((n_axial, self.n_channels))  # 液相轴向速度 (m/s)
        self.u_v = np.zeros((n_axial, self.n_channels))  # 气相轴向速度 (m/s)
        self.v_l = np.zeros((n_axial, self.n_channels))  # 液相横向速度 (m/s)
        self.v_v = np.zeros((n_axial, self.n_channels))  # 气相横向速度 (m/s)
        self.w_l = np.zeros((n_axial, self.n_channels))  # 液相z方向速度 (m/s)
        self.w_v = np.zeros((n_axial, self.n_channels))  # 气相z方向速度 (m/s)
        
        # 焓和温度场
        self.h_l = np.zeros((n_axial, self.n_channels))  # 液相焓 (J/kg)
        self.h_v = np.zeros((n_axial, self.n_channels))  # 气相焓 (J/kg)
        self.T_l = np.zeros((n_axial, self.n_channels))  # 液相温度 (K)
        self.T_v = np.zeros((n_axial, self.n_channels))  # 气相温度 (K)
        self.T_c = np.zeros((n_axial, self.n_channels))  # 包壳温度 (K)
        
        # 相变和传热参数
        self.gamma = np.zeros((n_axial, self.n_channels))  # 相变质量流量 (kg/m³·s)
        self.q_wl = np.zeros((n_axial, self.n_channels))  # 液相壁面热流密度 (W/m²)
        self.q_wv = np.zeros((n_axial, self.n_channels))  # 气相壁面热流密度 (W/m²)
        
        # === 瞬态计算参数 ===
        self.total_time = 5.0         # 总计算时间 (s)
        self.time_step = 0.001        # 时间步长 (s)
        self.current_time = 0.0       # 当前时间 (s)
        self.max_iterations = 1000    # 每个时间步内的最大迭代次数
        self.convergence_criteria = 1e-6  # 收敛判据
        
        # === 保存上一时间步的场值 ===
        self.prev_void_fraction = np.zeros((n_axial, self.n_channels))
        self.prev_liquid_velocity = np.zeros((n_axial, self.n_channels))
        self.prev_vapor_velocity = np.zeros((n_axial, self.n_channels))
        self.prev_pressure = np.zeros((n_axial, self.n_channels))
        self.prev_temperature = np.zeros((n_axial, self.n_channels))
        
        # === 初始化物理参数 ===
        self._initialize_physical_parameters()
        
        # === 结果保存参数 ===
        self.save_interval = 100  # 保存间隔步数
        self.results_folder = 'transient_results'  # 结果保存文件夹
        
        # 初始化所有场变量
        self.void_fraction = np.zeros((n_axial, self.n_channels))
        self.pressure = np.full((n_axial, self.n_channels), self.inlet_pressure)
        self.temperature = np.full((n_axial, self.n_channels), self.inlet_temp)
        self.liquid_velocity = np.zeros((n_axial, self.n_channels))
        self.vapor_velocity = np.zeros((n_axial, self.n_channels))
        self.T_l = np.full((n_axial, self.n_channels), self.inlet_temp)
        self.T_v = np.full((n_axial, self.n_channels), self.inlet_temp)
        self.T_c = np.full((n_axial, self.n_channels), self.inlet_temp)
        self.q_wl = np.zeros((n_axial, self.n_channels))
        self.q_wv = np.zeros((n_axial, self.n_channels))
        self.gamma = np.zeros((n_axial, self.n_channels))
        
        # 初始化历史数据字典
        self.history = {
            # 时间相关
            'time': [],
            'time_steps': [],
            
            # 场变量历史
            'void_fraction_history': [],
            'pressure_history': [],
            'temperature_history': [],
            'liquid_velocity_history': [],
            'vapor_velocity_history': [],
            'liquid_temp_history': [],
            'vapor_temp_history': [],
            'wall_temp_history': [],
            'heat_flux_liquid_history': [],
            'heat_flux_vapor_history': [],
            'phase_change_rate_history': [],
            
            # 统计量历史
            'max_void_fraction': [],
            'avg_pressure': [],
            'max_temperature': [],
            'max_liquid_velocity': [],
            'max_vapor_velocity': [],
            'residuals': [],
            'iterations': []
        }
        
        # 在初始化部分添加 w_prime_ij 数组
        self.w_prime_ij = np.zeros((self.n_channels, self.n_channels))  # 湍流混合系数矩阵
        
        # 初始化横向流动和混合系数矩阵
        self.w_ij = np.zeros((self.n_channels, self.n_channels))  # 横向流动系数矩阵
        
        # 计算通道直径
        self.channel_diameter = 2 * np.sqrt(self.channel_geometry[self.channel_types[1]]['flow_area'] / np.pi)
        
        # === 添加界面张力模型参数 ===
        self.surface_tension = 0.072  # 水的表面张力系数 (N/m)
        self.contact_angle = 45  # 接触角 (度)
        
        # === 添加湍流模型参数 (k-ε模型) ===
        self.k = np.zeros((n_axial, self.n_channels))  # 湍流动能
        self.epsilon = np.zeros((n_axial, self.n_channels))  # 湍流耗散率
        self.C_mu = 0.09  # k-ε模型常数
        self.C_1 = 1.44  # k-ε模型常数
        self.C_2 = 1.92  # k-ε模型常数
        self.sigma_k = 1.0  # k的湍流Schmidt数
        self.sigma_epsilon = 1.3  # ε的湍流Schmidt数
        
        # === 添加热耦合效应参数 ===
        self.T_fuel = np.zeros((n_axial, self.n_channels))  # 燃料温度场
        self.T_clad = np.zeros((n_axial, self.n_channels))  # 包壳温度场
        self.k_fuel = 3.0  # 燃料导热系数 (W/m·K)
        self.k_clad = 16.0  # 包壳导热系数 (W/m·K)
        self.gap_conductance = 10000.0  # 间隙导热系数 (W/m²·K)

        # 添加数值计算控制参数
        self.min_dt = 1e-8  # 最小时间步长
        self.max_dt = 1e-2  # 最大时间步长
        self.max_gamma = 1e3  # 最大相变率
        self.dynamic_tolerance = 1e-4  # 动态调整时间步长的误差阈值

    def initialize_field_variables(self):
        """
        初始化16子通道的场变量
        
        该��法初始���所有需要的物理��变量数组，包括:
        - 空泡率场
        - 压力场
        - 温度场
        - 速度场
        - 焓场
        - 传热相关场
        
        所有数组的大小为(n_axial, n_channels)
        """
        # 创建场变量数组
        self.void_fraction = np.zeros((self.n_axial, self.n_channels))
        self.pressure = np.zeros((self.n_axial, self.n_channels))
        self.temperature = np.zeros((self.n_axial, self.n_channels))
        self.liquid_velocity = np.zeros((self.n_axial, self.n_channels))
        self.vapor_velocity = np.zeros((self.n_axial, self.n_channels))
        self.T_l = np.zeros((self.n_axial, self.n_channels))
        self.T_v = np.zeros((self.n_axial, self.n_channels))
        self.T_c = np.zeros((self.n_axial, self.n_channels))
        self.h_l = np.zeros((self.n_axial, self.n_channels))
        self.h_v = np.zeros((self.n_axial, self.n_channels))
        self.q_wl = np.zeros((self.n_axial, self.n_channels))
        self.q_wv = np.zeros((self.n_axial, self.n_channels))
        self.gamma = np.zeros((self.n_axial, self.n_channels))

    def _initialize_physical_parameters(self):
        """
        初始化物理参数，考虑16子通道的特性
        
        该方法完成以下初始化工作�����
        1. 设置入口边界条件
        2. 初始化��向分布
        3. 为��同类型的子通道设置合适的初始值
        4. 考虑重力效应
        
        注意：
            - 角通道使用0.7的流量系数
            - 边通道使用0.9的流量系数
            - 中心通道使用1.0的流量系数
        """
        # 设置入口边界条件
        for j in range(self.n_channels):
            channel_type = self.channel_types[j+1]
            flow_area = self.channel_geometry[channel_type]['flow_area']
            
            # 入口空泡率和相速度
            self.void_fraction[0, j] = 0.0  # 入口无蒸汽
            self.alpha_l[0, j] = 1.0
            self.alpha_v[0, j] = 0.0
            
            # 入口压力和温度
            self.pressure[0, j] = self.inlet_pressure
            self.temperature[0, j] = self.inlet_temp
            self.T_l[0, j] = self.inlet_temp
            self.T_v[0, j] = self.inlet_temp
            self.T_c[0, j] = self.inlet_temp
            
            # 根据子通道类型计算入口速度
            rho_l = self.get_fluid_properties(self.inlet_pressure, self.inlet_temp, 'liquid')['density']
            # 根据子通道类型调整质量流量
            if channel_type == 'corner':
                mass_flux_factor = 0.7  # 角通道流量系数
            elif channel_type == 'edge':
                mass_flux_factor = 0.9  # 边通道流量系数
            else:  # center
                mass_flux_factor = 1.0  # 中心通道流量系数
                
            # 计算入口速度
            self.liquid_velocity[0, j] = (self.inlet_mass_flux * mass_flux_factor) / (rho_l * flow_area)
            self.vapor_velocity[0, j] = 0.0
            
            # 计算入口焓
            props_l = self.get_fluid_properties(self.inlet_pressure, self.inlet_temp, 'liquid')
            props_v = self.get_fluid_properties(self.inlet_pressure, self.inlet_temp, 'vapor')
            self.h_l[0, j] = props_l['enthalpy']
            self.h_v[0, j] = props_v['enthalpy']

        # 初始化轴向分布
        for i in range(1, self.n_axial):
            # 使用入口条件初始化整个流道
            self.void_fraction[i, :] = self.void_fraction[0, :]
            self.pressure[i, :] = self.pressure[0, :] - i * self.dz * 9.81 * 1000  # 考虑重力压降
            self.temperature[i, :] = self.temperature[0, :]
            self.liquid_velocity[i, :] = self.liquid_velocity[0, :]
            self.vapor_velocity[i, :] = self.vapor_velocity[0, :]
            self.T_l[i, :] = self.T_l[0, :]
            self.T_v[i, :] = self.T_v[0, :]
            self.T_c[i, :] = self.T_c[0, :]
            self.h_l[i, :] = self.h_l[0, :]
            self.h_v[i, :] = self.h_v[0, :]

    def calculate_initial_mass_flux_distribution(self):
        """
        计算初始质量流量分布
        
        根据子通道类型计算每个子通道的初始质量流量分布。
        考虑不同类型子通道的流动特性，使用不同的流量系数。
        
        返回：
            numpy.ndarray: 每个子通道的质量流量分布 (kg/m²·s)
        """
        total_flow_area = 0
        for j in range(self.n_channels):
            channel_type = self.channel_types[j+1]
            total_flow_area += self.channel_geometry[channel_type]['flow_area']
        
        # 计算平均质量流量
        avg_mass_flux = self.total_mass_flow / total_flow_area
        
        # 为每个子通道分配质量流量
        mass_flux_distribution = np.zeros(self.n_channels)
        for j in range(self.n_channels):
            channel_type = self.channel_types[j+1]
            if channel_type == 'corner':
                mass_flux_distribution[j] = avg_mass_flux * 0.7
            elif channel_type == 'edge':
                mass_flux_distribution[j] = avg_mass_flux * 0.9
            else:  # center
                mass_flux_distribution[j] = avg_mass_flux * 1.0
        
        return mass_flux_distribution

    def initialize_heat_flux(self):
        """
        初始化各子通道热流密度分布
        
        该方法：
        1. 考虑轴向功率分布（采用余弦分布）
        2. 考虑径向功率分布（不同类型子通道使用不同的功率因子）
        3. 初始时假设所有热量都传给液相
        
        功率因子：
        - 角通道：0.8
        - 边通道：0.9
        - 中心通道：1.0
        """
        for j in range(self.n_channels):
            channel_type = self.channel_types[j+1]
            heated_perimeter = self.channel_geometry[channel_type]['heated_perimeter']
            
            # 考虑轴向功率分布
            for i in range(self.n_axial):
                z = i * self.dz
                # 余弦分布的轴向功率因子
                axial_power_factor = np.cos(np.pi * z / (self.n_axial * self.dz))
                
                # 根据子通道类型设置径向功率因子
                if channel_type == 'corner':
                    radial_power_factor = 0.8
                elif channel_type == 'edge':
                    radial_power_factor = 0.9
                else:  # center
                    radial_power_factor = 1.0
                    
                # 计算局部热流密度
                q_base = self.average_heat_flux * axial_power_factor * radial_power_factor
                self.q_wl[i, j] = q_base * heated_perimeter
                self.q_wv[i, j] = 0.0  # 初始时假设所有热量都传给液相

    def get_fluid_properties(self, p, t, phase):
        """获取指定状态点的流体物性"""
        try:
            # 确保参数在CoolProp的有效范围内
            p = np.clip(p, 611.657, 22.064e6)  # 水的三相点到临界点压力范围
            t = np.clip(t, 273.16, 647.096)    # 水的三相点到临界点温度范围
            
            if phase == 'liquid':
                try:
                    return {
                        'density': PropsSI('D', 'P', p, 'T', t, 'Water'),
                        'viscosity': PropsSI('V', 'P', p, 'T', t, 'Water'),
                        'enthalpy': PropsSI('H', 'P', p, 'T', t, 'Water'),
                        'specific_heat': PropsSI('C', 'P', p, 'T', t, 'Water'),
                        'thermal_conductivity': PropsSI('L', 'P', p, 'T', t, 'Water')
                    }
                except Exception as e:
                    logging.error(f"液相物性计算失败 (P={p}, T={t}): {str(e)}")
                    return None
            else:  # vapor
                try:
                    return {
                        'density': PropsSI('D', 'P', p, 'Q', 1, 'Water'),
                        'viscosity': PropsSI('V', 'P', p, 'Q', 1, 'Water'),
                        'enthalpy': PropsSI('H', 'P', p, 'Q', 1, 'Water'),
                        'specific_heat': PropsSI('C', 'P', p, 'Q', 1, 'Water'),
                        'thermal_conductivity': PropsSI('L', 'P', p, 'Q', 1, 'Water')
                    }
                except Exception as e:
                    logging.error(f"气相物性计算失败 (P={p}, T={t}): {str(e)}")
                    return None
                
        except Exception as e:
            logging.error(f"物性计算失败 (P={p}, T={t}, phase={phase}): {str(e)}")
            return None

    def calculate_cross_flow(self, i, j, k):
        """
        计算子通道j和k之间的横向流动
        
        基于压力差和几何特征计算相邻子通道间的横向质量流量。
        考虑不同类型子通道的阻力特性。
        
        参数:
            i (int): 轴向位置索引
            j (int): 源子通道编号
            k (int): 目标子通道编号
            
        返回:
            float: 横向质量流量 (kg/s)，正值表示从j到k的流动
            
        注意:
            - 角通道的阻力系数为1.5
            - 边通道的阻力系数为1.2
            - 中心通道的阻力系数为1.0
        """
        try:
            # 检查子通道连接是否存在
            if k not in self.channel_connections[j]:
                return 0.0
            
            # 获取连接宽度
            gap_width = self.connection_gaps.get((j, k), 0.0)
            
            # 计算压力差
            dp = self.pressure[i, j-1] - self.pressure[i, k-1]
            
            # 获取子通道类型
            j_type = self.channel_types[j]
            k_type = self.channel_types[k]
            
            # 计算横向流动阻力系数
            if j_type == 'corner' or k_type == 'corner':
                resistance_factor = 1.5  # 角通道的阻力系数较大
            elif j_type == 'edge' or k_type == 'edge':
                resistance_factor = 1.2  # 边通道的阻力系数中等
            else:
                resistance_factor = 1.0  # 中心通道的阻力系数标准值
            
            # 计算横向流动系数
            s_jk = gap_width / (resistance_factor * self.dz)
            
            # 计算横向质量流量
            w_jk = s_jk * dp
            
            return w_jk
            
        except Exception as e:
            print(f"计算横向流动时出错: {str(e)}")
            return 0.0

    def calculate_mixing_coefficient(self, i, j, k):
        """
        计算子通道j和k之间的湍流混合系数
        
        基于流动特性和子通道类型计算湍流混合系数。
        考虑不同类型子通道的混合特性。
        
        参数：
            i (int): 轴向位置索引
            j (int): 第一个子通道编号
            k (int): 第二个子通道编号
            
        返回：
            float: 湍流混合系数 (无量纲)
            
        注意：
            - 角通道混合系数为基准值的0.7倍
            - 边通道混合系数为基准值的0.85倍
            - 中心通道混合系数为基准值
        """
        try:
            # 获取子通道类型
            j_type = self.channel_types[j]
            k_type = self.channel_types[k]
            
            # 获取平均流速
            v_j = self.liquid_velocity[i, j-1]
            v_k = self.liquid_velocity[i, k-1]
            v_avg = 0.5 * (v_j + v_k)
            
            # 基础混合系数
            beta_0 = 0.02
            
            # 根据子通道类型调整混合系数
            if j_type == 'corner' or k_type == 'corner':
                type_factor = 0.7  # 角通道混合较弱
            elif j_type == 'edge' or k_type == 'edge':
                type_factor = 0.85  # 边通道混合中等
            else:
                type_factor = 1.0  # 中心通道混合正常
            
            # 计算实际混合系数
            beta = beta_0 * type_factor * v_avg
            
            return beta
            
        except Exception as e:
            print(f"计算混合系数时出错: {str(e)}")
            return 0.0

    def calculate_turbulent_mixing(self, i, j, k, phi):
        """
        计算子通道j和k之间的湍流混合传输
        
        基于浓度梯度和混合系数计算相邻子通道间的湍流混合传输。
        
        参数:
            i (int): 轴向位置索引
            j, k (int): 相邻子通道编号
            phi (numpy.ndarray): 被混合的物理量（温度、焓等）
            
        返回:
            float: 混合传输量，正值表示从k到j的传输
            
        注意：
            - 混合方向由浓度梯度决定
            - 混合强度与连接宽度和浓度差值成正比
        """
        try:
            # 检查子通道连接
            if k not in self.channel_connections[j]:
                return 0.0
                
            # 获取连接宽度
            gap_width = self.connection_gaps.get((j, k), 0.0)
            
            # 计算混合系数
            beta = self.calculate_mixing_coefficient(i, j, k)
            
            # 计算物理量差值
            dphi = phi[i, j-1] - phi[i, k-1]
            
            # 计算混合流量
            w_prime = beta * gap_width * abs(dphi)
            
            # 确定混合方向
            if dphi > 0:
                return -w_prime  # 从j到k的混合
            else:
                return w_prime   # 从k到j的混合
                
        except Exception as e:
            print(f"计算湍流混合时出错: {str(e)}")
            return 0.0

    def update_cross_flow_terms(self, i):
        """
        更新所有子通道间的横向流动项
        """
        try:
            # 初始化横向流动数组
            self.w_ij = np.zeros((self.n_channels, self.n_channels))
            self.w_prime_ij = np.zeros((self.n_channels, self.n_channels))
            
            # 遍历所有子通道对
            for j in range(1, self.n_channels + 1):
                for k in self.channel_connections[j]:
                    if k > j:  # 避免重复计算
                        # 计算横向流动
                        w_jk = self.calculate_cross_flow(i, j, k)
                        self.w_ij[j-1, k-1] = w_jk
                        self.w_ij[k-1, j-1] = -w_jk
                        
                        # 计算湍流混合
                        # 温度混合
                        w_prime_t = self.calculate_turbulent_mixing(i, j, k, self.temperature)
                        # 焓混合
                        w_prime_h = self.calculate_turbulent_mixing(i, j, k, self.h_l)
                        
                        # 取平均值作为最终的湍流混合项
                        w_prime = 0.5 * (w_prime_t + w_prime_h)
                        self.w_prime_ij[j-1, k-1] = w_prime
                        self.w_prime_ij[k-1, j-1] = -w_prime
                        
        except Exception as e:
            print(f"更新横向流动项时出错: {str(e)}")

    def solve_mass_momentum(self):
        """求解质量和动量方程"""
        print("\n=== 开始求解质量和动量方程 ===")
        dt = self.dt  # 使用类的时间步长
        max_iterations = self.max_iterations
        tolerance = self.convergence_criteria
        
        for n in range(max_iterations):
            max_error = 0.0
            
            # 保存旧值用于计算误差
            old_void = np.copy(self.void_fraction)
            old_vl = np.copy(self.liquid_velocity)
            old_vv = np.copy(self.vapor_velocity)
            
            print(f"\n迭代步数: {n}")
            print(f"当前最大空泡率: {np.max(self.void_fraction):.4f}")
            print(f"当前最大液相速度: {np.max(self.liquid_velocity):.4f} m/s")
            print(f"当前最大气相速度: {np.max(self.vapor_velocity):.4f} m/s")
            
            # ���新横向流动项
            self.update_cross_flow_terms(n)
            
            for i in range(1, self.n_axial):
                for j in range(self.n_channels):
                    # 获���子通道类型和几何参数
                    channel_type = self.channel_types[j+1]
                    flow_area = self.channel_geometry[channel_type]['flow_area']
                    wetted_perimeter = self.channel_geometry[channel_type]['wetted_perimeter']
                    
                    print(f"\n处理子通道 {j+1} (类型: {channel_type})")
                    print(f"轴向位置: {i*self.dz:.3f} m")
                    print(f"流通面积: {flow_area:.6f} m²")
                    print(f"湿周长: {wetted_perimeter:.6f} m")
                    
                    # 获取流体属性
                    props_l = self.get_fluid_properties(self.pressure[i,j], self.temperature[i,j], 'liquid')
                    props_v = self.get_fluid_properties(self.pressure[i,j], self.temperature[i,j], 'vapor')
                    
                    if props_l is None or props_v is None:
                        print("警告: 无法获取流体物性")
                        continue
                        
                    rho_l = props_l['density']
                    rho_v = props_v['density']
                    
                    print(f"液相密度: {rho_l:.2f} kg/m³")
                    print(f"气相密度: {rho_v:.2f} kg/m³")
                    
                    # 计算横向流动的净贡献
                    net_cross_flow = 0.0
                    net_mixing = 0.0
                    for k in self.channel_connections[j+1]:
                        net_cross_flow += self.w_ij[j, k-1]
                        net_mixing += self.w_prime_ij[j, k-1]
                    
                    print(f"横向流动净贡献: {net_cross_flow:.4f} kg/s")
                    print(f"湍流混合净贡献: {net_mixing:.4f} kg/s")
                    
                    # 计算摩擦压降
                    f_l = self.calculate_friction_factor(self.liquid_velocity[i,j], 
                                                            props_l['viscosity'],
                                                        flow_area,
                                                        wetted_perimeter)
                                                    
                    f_v = self.calculate_friction_factor(self.vapor_velocity[i,j],
                                                            props_v['viscosity'],
                                                        flow_area,
                                                        wetted_perimeter)
                    
                    print(f"液相摩擦因子: {f_l:.4f}")
                    print(f"气相摩擦因子: {f_v:.4f}")
                    
                    # 使用原有的数值方法更新场变量
                    alpha = self.void_fraction[i-1,j]
                    vl = self.liquid_velocity[i-1,j]
                    vv = self.vapor_velocity[i-1,j]
                    
                    # 空泡率方程
                    new_alpha = alpha + dt * (
                        -vv * (alpha - self.void_fraction[i-1,j]) / self.dz
                        + net_cross_flow / (rho_v * flow_area)
                    )
                    new_alpha = np.clip(new_alpha, 0, 0.9999)
                    
                    print(f"更新后的空泡率: {new_alpha:.4f}")
                    
                    # 液相动量方程
                    new_vl = vl + dt * (
                        -vl * (vl - self.liquid_velocity[i-1,j]) / self.dz
                        - (1/rho_l) * (self.pressure[i,j] - self.pressure[i-1,j]) / self.dz
                        - f_l * vl * abs(vl) / (2 * self.channel_diameter)
                        - 9.81
                    )
                    
                    print(f"更新后的液相速度: {new_vl:.4f} m/s")
                    
                    # 气相动量方程
                    new_vv = vv + dt * (
                        -vv * (vv - self.vapor_velocity[i-1,j]) / self.dz
                        - (1/rho_v) * (self.pressure[i,j] - self.pressure[i-1,j]) / self.dz
                        - f_v * vv * abs(vv) / (2 * self.channel_diameter)
                        - 9.81
                    )
                    
                    print(f"更新后的气相速度: {new_vv:.4f} m/s")
                    
                    # 更新场变量
                    self.void_fraction[i,j] = new_alpha
                    self.liquid_velocity[i,j] = new_vl
                    self.vapor_velocity[i,j] = new_vv
                    
                    # 计算最大误差
                    error = max(
                        abs(new_alpha - old_void[i,j]),
                        abs(new_vl - old_vl[i,j]),
                        abs(new_vv - old_vv[i,j])
                    )
                    
                    max_error = max(max_error, error)
            
            print(f"\n当前迭代最大误差: {max_error:.6f}")
            
            # 检查收敛性并动态调整时间步长
            if max_error > self.dynamic_tolerance:
                dt *= 0.5
                print(f"误差过大，时间步长调整为 {dt:.2e}")
                if dt < self.min_dt:
                    print("错误：时间步长过小，计算终止")
                    break
            elif max_error < tolerance:
                print(f"\n计算在第{n}次迭代后收敛")
                break
            else:
                dt = min(dt * 1.1, self.max_dt)  # 适当增加时间步长但不超过最大值
            
            self.dt = dt
    
    def run_simulation(self):
        """运行完整的模拟计算"""
        print("开始数值计算...")
        
        time = 0
        step = 0
        
        while time < self.total_time:
            print(f"当前时间：{time:.3f} s")
            
            # 使用自适应时间步长解
            success = self.solve_two_fluid_equations()
            
            if not success:
                print(f"警告：在时间 {time:.3f} s 处计算不稳定")
                self.time_step *= 0.1
                if self.time_step < 1e-8:
                    print("错误：时间步长过小，计算终止")
                    break
                continue
            
            # 更新时间
            time += self.time_step
            step += 1
            
            # 保存结果
            if step % self.save_interval == 0:
                self.save_timestep_results(step)
                
        print("计算完成")
    
    def plot_results(self):
        """绘制计算结果的可视化图，增加错误检查和空值处理"""
        try:
            # 创建网格点
            x = np.arange(self.n_axial)
            y = np.arange(self.n_channels)
            X, Y = np.meshgrid(x, y)
            
            # 检查数据是否为空
            if np.all(self.void_fraction == 0) or np.all(self.pressure == 0):
                print("警告：检测到空数据，请确保计算已完成")
                return
                
            # 创建两个图形窗口
            fig1 = plt.figure(figsize=(20, 15))
            fig2 = plt.figure(figsize=(20, 15))
            
            # 第1个图形窗口：液相参数
            axes1 = []
            for i in range(6):
                ax = fig1.add_subplot(2, 3, i+1, projection='3d')
                axes1.append(ax)
                
            # 第二个图形窗口：气相参数
            axes2 = []
            for i in range(6):
                ax = fig2.add_subplot(2, 3, i+1, projection='3d')
                axes2.append(ax)
                
            # 绘制液相参数
            plots = [
                (self.alpha_l, '液相分数分布'),
                (self.u_l, '液相轴向速度分布 (m/s)'),
                (self.v_l, '液相横向速度分布 (m/s)'),
                (self.T_l, '液相温度分布 (K)'),
                (self.pressure, '压力分布 (Pa)'),
                (self.q_wl if hasattr(self, 'q_wl') else np.zeros_like(self.pressure), 
                 '液相壁面热流密度 (W/m²)')
            ]
            
            for ax, (data, title) in zip(axes1, plots):
                surf = ax.plot_surface(X, Y, data.T, cmap=cm.viridis)
                ax.set_title(title)
                ax.set_xlabel('轴向位置')
                ax.set_ylabel('子通道编号')
                fig1.colorbar(surf, ax=ax, shrink=0.5, aspect=5)
                
            # 绘制气相参数
            plots = [
                (self.alpha_v, '气相分数分布'),
                (self.u_v, '气相轴向速度分布 (m/s)'),
                (self.v_v, '气相横向速度分布 (m/s)'),
                (self.T_v, '气相温度分布 (K)'),
                (self.gamma if hasattr(self, 'gamma') else np.zeros_like(self.pressure), 
                 '相变率 (kg/m³·s)'),
                (self.q_wv if hasattr(self, 'q_wv') else np.zeros_like(self.pressure), 
                 '气相壁面热流密度 (W/m²)')
            ]
            
            for ax, (data, title) in zip(axes2, plots):
                surf = ax.plot_surface(X, Y, data.T, cmap=cm.viridis)
                ax.set_title(title)
                ax.set_xlabel('轴向位置')
                ax.set_ylabel('子通道编号')
                fig2.colorbar(surf, ax=ax, shrink=0.5, aspect=5)
                
            # 设置视角和网格
            for ax in axes1 + axes2:
                ax.view_init(elev=20, azim=45)
                ax.grid(True)
                
            plt.tight_layout()
            plt.show()
            
        except Exception as e:
            print(f"绘图时出错: {str(e)}")

    def solve_two_fluid_equations(self):
        """求解双流体模型方程组 - 添加自适应时间步长控制"""
        # 准备输入数据
        dx = dy = self.dz  # 假设网格等间距
        max_attempts = 5  # 最大重试次数
        
        for attempt in range(max_attempts):
            # 调用Numba优化的求解器
            results = solve_two_fluid_equations_numba(
                self.n_axial, self.n_channels,
                self.alpha_l, self.alpha_v,
                self.rho_l, self.rho_v,
                self.u_l, self.u_v,
                self.v_l, self.v_v,
                self.w_l, self.w_v,
                self.pressure, self.T_l, self.T_v, self.T_c,
                self.h_l, self.h_v,
                dx, dy, self.dz, self.time_step
            )
            
            # 解包结果
            if len(results) == 10:  # 正常返回
                (new_alpha_l, new_alpha_v, new_u_l, new_u_v,
                 new_v_l, new_v_v, new_h_l, new_h_v,
                 new_dt, needs_restart) = results
                
                if needs_restart:
                    print(f"时间��长调整为 {new_dt:.2e}，重新计算...")
                    self.time_step = new_dt
                    continue
                    
                # 更新时间步长
                self.time_step = new_dt
                
                # 更新场值
                self.alpha_l = new_alpha_l
                self.alpha_v = new_alpha_v
                self.u_l = new_u_l
                self.u_v = new_u_v
                self.v_l = new_v_l
                self.v_v = new_v_v
                self.h_l = new_h_l
                self.h_v = new_h_v
                
                print(f"计算成功，当前时间步长: {self.time_step:.2e}")
                return True
                
        print("警告：达到最大重试次数，计算可能不稳定")
        return False

    def calculate_phase_change(self, i, j):
        """计算相变质量流量"""
        try:
            T_c = self.T_c[i,j]  # 包壳温度
            T_sat = self.get_saturation_temperature(self.pressure[i,j])
            T_l = self.T_l[i,j]  # 液相温度
            
            # 过冷沸腾
            if T_l < T_sat and T_c > T_sat:
                h_l = self.calculate_subcooled_boiling_htc(T_c, T_sat, T_l, self.pressure[i,j])
                Q_l = np.pi * self.channel_diameter * h_l * (T_c - T_l)
                gamma_l_sc = Q_l / self.get_latent_heat(self.pressure[i,j])
                gamma_l_sc = np.clip(gamma_l_sc, -self.max_gamma, 0)  # 限制相变率范围
                return gamma_l_sc
                
            # 饱和核态沸腾
            elif T_l >= T_sat and T_c > T_sat:
                h_l = self.calculate_nucleate_boiling_htc(T_c, T_sat, self.pressure[i,j])
                Q_l = np.pi * self.channel_diameter * h_l * (T_c - T_sat)
                gamma_l_sh = Q_l / self.get_latent_heat(self.pressure[i,j])
                gamma_l_sh = np.clip(gamma_l_sh, 0, self.max_gamma)  # 限制相变率范围
                return gamma_l_sh
                
            return 0.0
            
        except Exception as e:
            print(f"计算相变率时出错: {str(e)}")
            return 0.0

    def calculate_turbulent_mixing(self, i, j, phase):
        """计算湍流混合和空泡漂移质量流量"""
        C_e = 0.01  # 涡流扩散率
        L_tur = self.channel_diameter  # 湍流混合长度
        
        mixing_mass = 0.0
        for k in self.channel_connections[j+1]:
            L_gap = self.get_gap_width(j+1, k)  # 获取间隙宽度
            
            if phase == 'liquid':
                rho_diff = self.get_fluid_properties(self.pressure[i,j], 
                    self.temperature[i,j], 'liquid')['density']
                alpha_diff = self.alpha_l[i,j] - self.alpha_l[i,k-1]
                alpha_eq_diff = self.calculate_equilibrium_void_fraction(i,j) - \
                               self.calculate_equilibrium_void_fraction(i,k-1)
            else:
                rho_diff = self.get_fluid_properties(self.pressure[i,j],
                    self.temperature[i,j], 'vapor')['density']
                alpha_diff = self.alpha_v[i,j] - self.alpha_v[i,k-1]
                alpha_eq_diff = (1-self.calculate_equilibrium_void_fraction(i,j)) - \
                               (1-self.calculate_equilibrium_void_fraction(i,k-1))
                
            mixing_mass += C_e * L_gap * L_tur * rho_diff * (alpha_diff - alpha_eq_diff)
            
        return mixing_mass

    def calculate_wall_friction(self, i, j, phase, direction='axial'):
        """计算壁面摩擦力
        
        参数:
        i, j: 网格位置
        phase: 'liquid'或'vapor'
        direction: 'axial'或'lateral'
        """
        if phase == 'liquid':
            rho = self.get_fluid_properties(self.pressure[i,j], self.T_l[i,j], 'liquid')['density']
            mu = self.get_fluid_properties(self.pressure[i,j], self.T_l[i,j], 'liquid')['viscosity']
            if direction == 'axial':
                v = self.u_l[i,j]
            else:
                v = self.v_l[i,j]
        else:
            rho = self.get_fluid_properties(self.pressure[i,j], self.T_v[i,j], 'vapor')['density']
            mu = self.get_fluid_properties(self.pressure[i,j], self.T_v[i,j], 'vapor')['viscosity']
            if direction == 'axial':
                v = self.u_v[i,j]
            else:
                v = self.v_v[i,j]
            
        # 计算Reynolds数
        Re = rho * abs(v) * self.channel_diameter / mu
        
        # 计算摩擦系数
        if Re < 2300:  # 层流
            f = 64 / Re
        else:  # 湍流
            f = 0.316 * Re**(-0.25)
            
        # 计算壁面摩擦力
        tau_w = 0.5 * f * rho * v * abs(v)
        
        # 考虑两相流修正系数
        if phase == 'liquid':
            phi_l = (1 - self.alpha_l[i,j])**(-1.75)  # 两相流修正系数
            tau_w *= phi_l
            
        return tau_w * np.pi * self.channel_diameter / self.channel_area

    def calculate_interfacial_force(self, i, j, direction='axial'):
        """计算界面剪切力 (eq.19)"""
        rho_l = self.get_fluid_properties(self.pressure[i,j], self.T_l[i,j], 'liquid')['density']
        
        if direction == 'axial':
            v_rel = self.u_v[i,j] - self.u_l[i,j]  # 相对速度
        else:
            v_rel = self.v_v[i,j] - self.v_l[i,j]
            
        # 计算界面面积密度
        a_i = 6 * self.alpha_v[i,j] / self.channel_diameter  # 假设球形气泡
        
        # 计算界面阻力系数
        C_i = 0.44  # 球形气泡的阻力系数
        
        # 计算界面剪切力
        tau_i = 0.5 * C_i * rho_l * v_rel * abs(v_rel) * a_i
        
        return tau_i

    def calculate_wall_heat_flux(self, i, j, phase):
        """计算壁面热流密度 (eq.6)"""
        if phase == 'liquid':
            # 根据沸腾状态确定换热系数
            T_sat = self.get_saturation_temperature(self.pressure[i,j])
            
            if self.T_l[i,j] < T_sat and self.T_c[i,j] > T_sat:
                # 过冷核态沸腾
                h = self.calculate_subcooled_boiling_htc(
                    self.T_c[i,j], T_sat, self.T_l[i,j], self.pressure[i,j])
                
                
            elif self.T_l[i,j] >= T_sat and self.T_c[i,j] > T_sat:
                # 饱和核态沸腾
                h = self.calculate_nucleate_boiling_htc(
                    self.T_c[i,j], T_sat, self.pressure[i,j])
                
                
            else:
                # 单相对流换热
                h = self.calculate_single_phase_htc(i, j, 'liquid')
                
            q = h * (self.T_c[i,j] - self.T_l[i,j])
        else:
            # 气相单相对流换热
            h = self.calculate_single_phase_htc(i, j, 'vapor')
            q = h * (self.T_c[i,j] - self.T_v[i,j])
            
        return q * np.pi * self.channel_diameter

    def calculate_subcooled_boiling_htc(self, T_w, T_sat, T_l, P):
        """计算过冷核态沸腾换热系数 (eq.8)"""
        F_gam = 0.5  # 气泡冷凝系数
        h = (1 - F_gam) * np.exp(2*P/1260 * (T_w-T_sat)/(2722*(T_w-T_l)))
        return h

    def calculate_nucleate_boiling_htc(self, T_w, T_sat, P):
        """计算饱和核态沸腾换热系数 (eq.9)"""
        h = np.exp(2*P/1260 * (T_w-T_sat)/2722)
        return h

    def calculate_radial_heat_conduction(self, i, j):
        """计算径向热传导项"""
        # 简化处理,假设均匀网格
        dr = self.channel_diameter / 10  # 径向网格间距
        r = np.linspace(dr/2, self.channel_diameter/2, 10)
        
        # 计算径向温度梯度
        d2T_dr2 = np.zeros_like(r)
        for k in range(1, len(r)-1):
            d2T_dr2[k] = (1/r[k]) * (
                (r[k+1]*self.T_c[i,j] - r[k]*self.T_c[i,j])/dr -
                (r[k]*self.T_c[i,j] - r[k-1]*self.T_c[i,j])/dr)
        
        return np.mean(d2T_dr2)  # 返回平均值

    def calculate_fuel_heat_generation(self, i, j):
        """计算燃料热生成"""
        # 简化处理,假设轴向余弦分布
        q_max = 1e6  # 最大线功率密度 (W/m)
        z = i * self.dz
        H = self.channel_length
        q = q_max * np.cos(np.pi * z / H)
        return q / (np.pi * self.channel_diameter)

    def run_transient_simulation(self):
        """运行瞬态计算"""
        print("开始瞬态计算...")
        
        n_time_steps = int(self.total_time / self.time_step)
        
        for time_step in range(n_time_steps):
            # ���存上一时间步的值
            self.save_previous_step_values()
            
            # 在当前时间步内进行迭代求解
            converged, n_iter, final_residual = self.solve_current_time_step()
            
            # 记录残差和迭代次数
            self.history['residuals'].append(final_residual)
            self.history['iterations'].append(n_iter)
            
            # 保存结果
            self.save_timestep_results(time_step)
            
            if not converged:
                print(f"警告：时间步 {time_step + 1} 未达到收敛")
                
            # 更新当前时间
            self.current_time += self.time_step
            
            # 打印计算进度
            if (time_step + 1) % 100 == 0:
                print(f"计算进度: {(time_step + 1)/n_time_steps*100:.1f}% "
                      f"(当前时间: {self.current_time:.3f} s)")
        
        print("瞬态计算完成")
        print(f"结果已保存在 '{self.results_folder}' 文件夹中")
        
    def save_previous_step_values(self):
        """保存上一时间步的场值"""
        self.prev_void_fraction = np.copy(self.void_fraction)
        self.prev_liquid_velocity = np.copy(self.liquid_velocity)
        self.prev_vapor_velocity = np.copy(self.vapor_velocity)
        self.prev_pressure = np.copy(self.pressure)
        self.prev_temperature = np.copy(self.temperature)

    def solve_current_time_step(self):
        """在当前时间步内进行迭代求解"""
        residuals = []
        
        for iteration in range(self.max_iterations):
            # 保存当前迭代步的值用于收敛性检查
            old_void = np.copy(self.void_fraction)
            old_vl = np.copy(self.liquid_velocity)
            old_vv = np.copy(self.vapor_velocity)
            old_p = np.copy(self.pressure)
            old_t = np.copy(self.temperature)
            
            # 求解质量和动量方程
            self.solve_mass_momentum_equations()
            
            # 求解能量方程
            self.solve_energy_equations()
            
            # 检查收敛性
            max_error = self.check_convergence(
                old_void, old_vl, old_vv, old_p, old_t)
            
            residuals.append(max_error)
            
            if max_error < self.convergence_criteria:
                return True, iteration + 1, max_error
                
        return False, self.max_iterations, residuals[-1]
    
    def save_timestep_results(self, time_step):
        """保存当前时间步的结果"""
        try:
            if time_step % self.save_interval == 0:
                # 1. 更新时间相关数据
                self.history['time'].append(self.current_time)
                self.history['time_steps'].append(time_step)
                
                # 2. 更新场变量历史
                self.history['void_fraction_history'].append(np.copy(self.void_fraction))
                self.history['pressure_history'].append(np.copy(self.pressure))
                self.history['temperature_history'].append(np.copy(self.temperature))
                self.history['liquid_velocity_history'].append(np.copy(self.liquid_velocity))
                self.history['vapor_velocity_history'].append(np.copy(self.vapor_velocity))
                self.history['liquid_temp_history'].append(np.copy(self.T_l))
                self.history['vapor_temp_history'].append(np.copy(self.T_v))
                self.history['wall_temp_history'].append(np.copy(self.T_c))
                self.history['heat_flux_liquid_history'].append(np.copy(self.q_wl))
                self.history['heat_flux_vapor_history'].append(np.copy(self.q_wv))
                self.history['phase_change_rate_history'].append(np.copy(self.gamma))
                
                # 3. 更新统计量历史
                self.history['max_void_fraction'].append(np.max(self.void_fraction))
                self.history['avg_pressure'].append(np.mean(self.pressure))
                self.history['max_temperature'].append(np.max(self.temperature))
                self.history['max_liquid_velocity'].append(np.max(self.liquid_velocity))
                self.history['max_vapor_velocity'].append(np.max(self.vapor_velocity))
                
                # 4. 保存CSV文件
                self.save_to_csv(time_step)
                
                # 5. 定期保存完整状态（每100个时间步）
                if time_step % (self.save_interval * 100) == 0:
                    self.save_full_state(time_step)
                    
        except Exception as e:
            print(f"保存时间步 {time_step} 的结果时出错: {str(e)}")
            print(f"错误类型: {type(e).__name__}")
            print(f"错误位置: {traceback.format_exc()}")
            
    def save_to_csv(self, time_step):
        """保存当前时间步的数据到CSV文件"""
        try:
            filename = os.path.join(self.results_folder, f"timestep_{time_step:06d}.csv")
            
            # 使用UTF-8编码保存
            with open(filename, 'w', newline='', encoding='utf-8') as csvfile:
                writer = csv.writer(csvfile)
                
                # 1. 写入时间步信息
                writer.writerow(['时间步信息'])
                writer.writerow(['当前时间 (s)', f"{self.current_time:.6f}"])
                writer.writerow(['时间步编号', time_step])
                writer.writerow([])
                
                # 2. 写入全局参数
                writer.writerow(['全局参数'])
                writer.writerow(['平均压力 (Pa)', f"{np.mean(self.pressure):.2f}"])
                writer.writerow(['最大温度 (K)', f"{np.max(self.temperature):.2f}"])
                writer.writerow(['最大空泡率', f"{np.max(self.void_fraction):.4f}"])
                writer.writerow([])
                
                # 3. 写入详细数据
                headers = [
                    '轴向位置 (m)',
                    '子通道编号',
                    '空泡率',
                    '压力 (Pa)',
                    '液相温度 (K)',
                    '气相温度 (K)',
                    '壁面温度 (K)',
                    '液相速度 (m/s)',
                    '气相速度 (m/s)',
                    '液相壁面热流密度 (W/m²)',
                    '气相壁面热流密度 (W/m²)',
                    '相变率 (kg/m³·s)'
                ]
                writer.writerow(headers)
                
                # 写入每个网格点的数据
                for i in range(self.n_axial):
                    for j in range(self.n_channels):
                        row_data = [
                            f"{i * self.dz:.6f}",
                            j + 1,
                            f"{self.void_fraction[i, j]:.6f}",
                            f"{self.pressure[i, j]:.2f}",
                            f"{self.T_l[i, j]:.2f}",
                            f"{self.T_v[i, j]:.2f}",
                            f"{self.T_c[i, j]:.2f}",
                            f"{self.liquid_velocity[i, j]:.4f}",
                            f"{self.vapor_velocity[i, j]:.4f}",
                            f"{self.q_wl[i, j]:.2f}",
                            f"{self.q_wv[i, j]:.2f}",
                            f"{self.gamma[i, j]:.6f}"
                        ]
                        writer.writerow(row_data)
                        
            print(f"结果已保存到: {filename}")
            
        except Exception as e:
            print(f"保存CSV文件时出错: {str(e)}")
            
    def save_full_state(self, time_step):
        """保存完整的计算状态"""
        try:
            filename = os.path.join(self.results_folder, f"full_state_{time_step:06d}.npz")
            np.savez_compressed(
                filename,
                time=self.current_time,
                time_step=time_step,
                void_fraction=self.void_fraction,
                pressure=self.pressure,
                temperature=self.temperature,
                liquid_velocity=self.liquid_velocity,
                vapor_velocity=self.vapor_velocity,
                T_l=self.T_l,
                T_v=self.T_v,
                T_c=self.T_c,
                q_wl=self.q_wl,
                q_wv=self.q_wv,
                gamma=self.gamma
            )
            print(f"完整状态已保存到: {filename}")
            
        except Exception as e:
            print(f"保存完整状态时出错: {str(e)}")

    def plot_transient_results(self):
        """绘制瞬态计算结果"""
        try:
            # 创建一个包含4个子图的图形
            fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12)) # 移除多余的括号
            
            # 1. 残差历史
            ax1.semilogy(self.history['residuals'])
            ax1.set_title('残差历史')
            ax1.set_xlabel('时间步')
            ax1.set_ylabel('残差')
            ax1.grid(True)
            
            # 2. 最大空泡率随时间的变化
            ax2.plot(self.history['time'], self.history['max_void_fraction'])
            ax2.set_title('最大空泡率随时间的变化')
            ax2.set_xlabel('时间 (s)')
            ax2.set_ylabel('最大空泡率')
            ax2.grid(True)
            
            # 3. 平均压力随时间的变化
            ax3.plot(self.history['time'], self.history['avg_pressure'])
            ax3.set_title('平均压力随时间的变化')
            ax3.set_xlabel('时间 (s)')
            ax3.set_ylabel('压力 (Pa)')
            ax3.grid(True)
            
            # 4. 最大温度随时间的变化
            ax4.plot(self.history['time'], self.history['max_temperature'])
            ax4.set_title('最大温度随时间的变化')
            ax4.set_xlabel('时间 (s)')
            ax4.set_ylabel('温度 (K)')
            ax4.grid(True)
            
            plt.tight_layout()
            plt.show()
            
            # 创建动画展示场分布随时间的变化
            self.animate_field_evolution()
            
        except Exception as e:
            print(f"绘制瞬态结果时出错: {str(e)}")

    def animate_field_evolution(self):
        """创建场分布随时间变化的动画并保存"""
        try:
            # 创建动画图形
            fig = plt.figure(figsize=(15, 5))
            
            # 创建三个子图：空泡率、温度和压力
            ax1 = fig.add_subplot(131, projection='3d')
            ax2 = fig.add_subplot(132, projection='3d')
            ax3 = fig.add_subplot(133, projection='3d')
            
            x = np.arange(self.n_axial)
            y = np.arange(self.n_channels)
            X, Y = np.meshgrid(x, y)
            
            # 确保有足够的历史数据
            n_frames = min(len(self.history['time']),
                      len(self.history['void_fraction_history']),
                      len(self.history['temperature_history']),
                      len(self.history['pressure_history']))
            
            if n_frames == 0:
                print("警告：没有可用的历史数据来创建动画")
                return
                
            def update(frame):
                # 清除当前帧
                ax1.clear()
                ax2.clear()
                ax3.clear()
                
                # 确保frame不会超出数组范围
                frame = min(frame, n_frames - 1)
                
                try:
                    # 绘制空泡率分布
                    surf1 = ax1.plot_surface(X, Y, 
                        self.history['void_fraction_history'][frame].T,
                        cmap=cm.viridis)
                    ax1.set_title(f'空泡率分布\n时间: {self.history["time"][frame]:.3f} s')
                    
                    # 绘制温度分布
                    surf2 = ax2.plot_surface(X, Y,
                        self.history['temperature_history'][frame].T,
                        cmap=cm.viridis)
                    ax2.set_title('温度分布 (K)')
                    
                    # 绘制压力分布
                    surf3 = ax3.plot_surface(X, Y,
                        self.history['pressure_history'][frame].T,
                        cmap=cm.viridis)
                    ax3.set_title('压力分布 (Pa)')
                    
                    # 设置视角
                    for ax in [ax1, ax2, ax3]:
                        ax.view_init(elev=20, azim=45)
                        ax.grid(True)
                        ax.set_xlabel('轴向位置')
                        ax.set_ylabel('子通道编号')
                        
                except IndexError as e:
                    print(f"警告：动画帧 {frame} 超出数据范围")
                    return
            
            # 创建动画
            anim = FuncAnimation(fig, update, 
                                frames=n_frames,
                                interval=200)
            
            # 保存动画
            try:
                anim.save(f'{self.results_folder}/evolution.mp4', 
                         writer='ffmpeg', fps=10)
                print(f"动画已保存到 {self.results_folder}/evolution.mp4")
            except Exception as e:
                print(f"保存动画失败: {str(e)}")
                print("尝试使用其他格式保存...")
                try:
                    anim.save(f'{self.results_folder}/evolution.gif', 
                             writer='pillow', fps=10)
                    print(f"动画已保存为GIF格式")
                except Exception as e:
                    print(f"保存GIF也失败: {str(e)}")
            
            plt.tight_layout()
            plt.show()
            
        except Exception as e:
            print(f"创建动画时出错: {str(e)}")

    def create_results_folder(self):
        """创建结果保存文件夹"""
        if not os.path.exists(self.results_folder):
            os.makedirs(self.results_folder)
            

    def solve_mass_momentum_equations(self):
        """求解质量和动量方程 - Numba优化版本"""
        # 使用Numba优化的网格点计算
        new_void, new_vl, new_vv = compute_grid_points_numba(
            self.n_axial, self.n_channels,
            self.void_fraction, self.prev_void_fraction,
            self.vapor_velocity, self.liquid_velocity,
            self.prev_liquid_velocity, self.prev_vapor_velocity,
            self.pressure, self.temperature,
            self.time_step, self.dz
        )
        
        # 更新结果
        self.void_fraction = new_void
        self.liquid_velocity = new_vl
        self.vapor_velocity = new_vv
        
    def solve_energy_equations(self):
        """求解能量方程，适应16子通道结构"""
        try:
            dt = 1e-4  # 时间步长
            
            for i in range(1, self.n_axial):
                for j in range(self.n_channels):
                    # 获取流体物性
                    props_l = self.get_fluid_properties(self.pressure[i,j], self.T_l[i,j], 'liquid')
                    props_v = self.get_fluid_properties(self.pressure[i,j], self.T_v[i,j], 'vapor')
                    
                    if props_l is None or props_v is None:
                        logging.warning(f"子通道 {j+1} 在位置 {i} 处无法获取流体物性")
                        continue
                    
                    # 获取饱和焓
                    try:
                        h_sat_l = PropsSI('H', 'P', self.pressure[i,j], 'Q', 0, 'Water')
                        h_sat_v = PropsSI('H', 'P', self.pressure[i,j], 'Q', 1, 'Water')
                    except Exception as e:
                        logging.error(f"计算饱和焓失败: {str(e)}")
                        continue
                    
                    # 计算相变率
                    gamma = self.calculate_phase_change_rate(i, j)
                    
                    # 计算液相焓变化
                    dh_l = dt * (
                        -self.liquid_velocity[i,j] * (self.h_l[i,j] - self.h_l[i-1,j]) / self.dz +
                    self.q_wl[i,j] / (props_l['density'] * (1 - self.void_fraction[i,j])) +
                    -gamma * self.h_l[i,j] / props_l['density']
                    )
                    # 计算气相焓变化
                    dh_v = dt * (
                        -self.vapor_velocity[i,j] * (self.h_v[i,j] - self.h_v[i-1,j]) / self.dz +
                    self.q_wv[i,j] / (props_v['density'] * self.void_fraction[i,j]) +
                    gamma * self.h_v[i,j] / props_v['density']
                    )
                    # 更新焓值并限制在物理范围内
                    self.h_l[i,j] = np.clip(self.h_l[i,j] + dh_l, props_l['enthalpy'], h_sat_v)
                    self.h_v[i,j] = np.clip(self.h_v[i,j] + dh_v, h_sat_l, h_sat_v * 1.2)
                    
                    # 更新温度
                    try:
                        self.T_l[i,j] = PropsSI('T', 'P', self.pressure[i,j], 'H', self.h_l[i,j], 'Water')
                        self.T_v[i,j] = PropsSI('T', 'P', self.pressure[i,j], 'H', self.h_v[i,j], 'Water')
                    except Exception as e:
                        logging.warning(f"温度计算失败，使用上一时间步的值")
                        self.T_l[i,j] = self.T_l[i-1,j]
                        self.T_v[i,j] = self.T_v[i-1,j]
                    
                    # 限制温度在物理范围内
                    self.T_l[i,j] = np.clip(self.T_l[i,j], 273.15, 647.0)
                    self.T_v[i,j] = np.clip(self.T_v[i,j], 273.15, 647.0)
                    
        except Exception as e:
            logging.error(f"求解能量方程时出错: {str(e)}")
            logging.error(f"错误位置: {traceback.format_exc()}")
    
    def calculate_phase_change_rate(self, i, j):
        """
        计算相变率，考虑多种相变机制
        
        包括:
        1. 过冷沸腾
        2. 饱和沸腾
        3. 凝结
        """
        try:
            # 获取局部条件
            P = self.pressure[i,j]
            T_l = self.T_l[i,j]
            T_v = self.T_v[i,j]
            T_w = self.T_c[i,j]
            
            # 使用CoolProp获取精确的热力学参数
            T_sat = PropsSI('T', 'P', P, 'Q', 0, 'Water')
            h_fg = PropsSI('H', 'P', P, 'Q', 1, 'Water') - PropsSI('H', 'P', P, 'Q', 0, 'Water')
            
            # 获取流体物性
            props_l = self.get_fluid_properties(P, T_l, 'liquid')
            props_v = self.get_fluid_properties(P, T_v, 'vapor')
            
            if props_l is None or props_v is None:
                return 0.0
                
            # 1. 过冷沸腾
            if T_l < T_sat and T_w > T_sat:
                # 使用改进的过冷沸腾模型
                q_sc = self.calculate_subcooled_boiling_heat_flux(T_w, T_sat, T_l, P)
                gamma_sc = q_sc / h_fg
                
            # 2. 饱和沸腾
            elif T_l >= T_sat and T_w > T_sat:
                # 使用改进的核态沸腾模型
                q_nb = self.calculate_nucleate_boiling_heat_flux(T_w, T_sat, P)
                gamma_nb = q_nb / h_fg
                
            # 3. 凝结
            elif T_v > T_sat:
                # 使用改进的凝结模型
                q_cond = self.calculate_condensation_heat_flux(T_v, T_sat, P)
                gamma_cond = -q_cond / h_fg
                
            else:
                return 0.0
                
            # 考虑界面面积密度的影响
            a_i = 6 * self.void_fraction[i,j] / 0.001  # 假设气泡直径为1mm
            
            # 计算总的相变率
            if T_l < T_sat and T_w > T_sat:
                gamma = gamma_sc * a_i
            elif T_l >= T_sat and T_w > T_sat:
                gamma = gamma_nb * a_i
            else:
                gamma = gamma_cond * a_i
                
            return gamma
            
        except Exception as e:
            print(f"计算相变率时出错: {str(e)}")
            return 0.0

    def calculate_subcooled_boiling_heat_flux(self, T_w, T_sat, T_l, P):
        """改进的过冷沸腾换热模型"""
        try:
            # 获取流体物性
            props_l = self.get_fluid_properties(P, T_l, 'liquid')
            if props_l is None:
                return 0.0
                
            # 计算对流换热系数
            h_c = self.calculate_single_phase_htc(T_w, P, 'liquid')
            
            # 计算沸腾抑制因子
            S = max(0, min(1, (T_w - T_sat)/(T_sat - T_l)))
            
            # 计算核态沸腾换热系数
            h_nb = self.calculate_nucleate_boiling_htc(T_w, T_sat, P)
            
            # 组合换热系数
            q_sc = h_c * (T_w - T_l) + S * h_nb * (T_w - T_sat)
            
            return q_sc
            
        except Exception as e:
            print(f"计算过冷沸腾热流密度时出错: {str(e)}")
            return 0.0

    def calculate_nucleate_boiling_heat_flux(self, T_w, T_sat, P):
        """改进的核态沸腾换热模型"""
        try:
            # 使用改进的Rohsenow相关式
            props_l = self.get_fluid_properties(P, T_sat, 'liquid')
            props_v = self.get_fluid_properties(P, T_sat, 'vapor')
            
            if props_l is None or props_v is None:
                return 0.0
                
            h_fg = PropsSI('H', 'P', P, 'Q', 1, 'Water') - PropsSI('H', 'P', P, 'Q', 0, 'Water')
            sigma = 0.072  # 水的表面张力系数
            g = 9.81
            
            # 计算无量纲参数
            Ja = props_l['specific_heat'] * (T_w - T_sat) / h_fg
            Ra = (g * (props_l['density'] - props_v['density']) * 
                  props_l['density']**2 * h_fg * self.channel_diameter**3 /
                  (props_l['viscosity'] * props_l['thermal_conductivity'] * (T_w - T_sat)))
                  
            # 计算核态沸腾热流密度
            q_nb = props_l['thermal_conductivity'] * (T_w - T_sat) / self.channel_diameter * \
                   (Ja * Ra)**0.33
                   
            return q_nb
            
        except Exception as e:
            print(f"计算核态沸腾热流密度时出错: {str(e)}")
            return 0.0

    def calculate_condensation_heat_flux(self, T_v, T_sat, P):
        """改进的凝结换热模型"""
        try:
            # 使用Nusselt凝结理论
            props_l = self.get_fluid_properties(P, T_sat, 'liquid')
            props_v = self.get_fluid_properties(P, T_sat, 'vapor')
            
            if props_l is None or props_v is None:
                return 0.0
                
            h_fg = PropsSI('H', 'P', P, 'Q', 1, 'Water') - PropsSI('H', 'P', P, 'Q', 0, 'Water')
            g = 9.81
            
            # ��算凝结膜厚度
            delta = (4 * props_l['viscosity'] * props_l['thermal_conductivity'] * 
                    (T_v - T_sat) * self.channel_length /
                    (g * props_l['density'] * (props_l['density'] - props_v['density']) * h_fg))**0.25
                    
            # 计算凝结热流密度
            q_cond = props_l['thermal_conductivity'] * (T_v - T_sat) / delta
            
            return q_cond
            
        except Exception as e:
            print(f"计算凝结热流密度时出错: {str(e)}")
            return 0.0

    def calculate_turbulent_mixing(self, i, j, phase):
        """计算湍流混合和空泡漂移质量流量"""
        C_e = 0.01  # 涡流扩散率
        L_tur = self.channel_diameter  # 湍流混合长度
        
        mixing_mass = 0.0
        for k in self.channel_connections[j+1]:
            L_gap = self.get_gap_width(j+1, k)  # 获取间隙宽度
            
            if phase == 'liquid':
                rho_diff = self.get_fluid_properties(self.pressure[i,j], 
                    self.temperature[i,j], 'liquid')['density']
                alpha_diff = self.alpha_l[i,j] - self.alpha_l[i,k-1]
                alpha_eq_diff = self.calculate_equilibrium_void_fraction(i,j) - \
                               self.calculate_equilibrium_void_fraction(i,k-1)
            else:
                rho_diff = self.get_fluid_properties(self.pressure[i,j],
                    self.temperature[i,j], 'vapor')['density']
                alpha_diff = self.alpha_v[i,j] - self.alpha_v[i,k-1]
                alpha_eq_diff = (1-self.calculate_equilibrium_void_fraction(i,j)) - \
                               (1-self.calculate_equilibrium_void_fraction(i,k-1))
                
            mixing_mass += C_e * L_gap * L_tur * rho_diff * (alpha_diff - alpha_eq_diff)
            
        return mixing_mass

    def calculate_wall_friction(self, i, j, phase, direction='axial'):
        """计算壁面摩擦力"""
        print(f"\n=== 计算壁面摩擦力 (i={i}, j={j}, phase={phase}, direction={direction}) ===")
        if phase == 'liquid':
            props = self.get_fluid_properties(self.pressure[i,j], self.T_l[i,j], 'liquid')
            if direction == 'axial':
                velocity = self.liquid_velocity[i,j]
            else:
                velocity = self.v_l[i,j]
        else:
            props = self.get_fluid_properties(self.pressure[i,j], self.T_v[i,j], 'vapor')
            if direction == 'axial':
                velocity = self.vapor_velocity[i,j]
            else:
                velocity = self.v_v[i,j]
                
        if props is None:
            print("警告: 无法获取流体物性")
            return 0.0
            
        print(f"流体速度: {velocity:.4f} m/s")
        print(f"流体密度: {props['density']:.2f} kg/m³")
        print(f"流体粘度: {props['viscosity']:.6f} Pa·s")
        
        friction = calculate_wall_friction_numba(
            velocity, self.channel_diameter,
            props['density'], props['viscosity']
        )
        
        print(f"计算得到的壁面摩擦力: {friction:.4f} N/m²")
        return friction

    def calculate_residual(self):
        """计算残差 - Numba优化版本"""
        return calculate_residual_numba(
            self.n_axial, self.n_channels,
            self.liquid_velocity, self.liquid_enthalpy,
            self.dz
        )

    def get_latent_heat(self, pressure):
        """获取给定压力下的汽化潜热"""
        try:
            h_f = PropsSI('H', 'P', pressure, 'Q', 0, 'Water')  # 饱和液焓
            h_g = PropsSI('H', 'P', pressure, 'Q', 1, 'Water')  # 饱和气焓
            return h_g - h_f
        except:
            return 2.257e6  # 返回1个大气压下的近似值

    def get_saturation_temperature(self, pressure):
        """
        计算给定压力下的饱和温度
        
        参数：
        pressure: float, 压力 (Pa)
        
        返回：
        float: 饱和温度 (K)
        """
        try:
            return PropsSI('T', 'P', pressure, 'Q', 0, 'Water')
        except Exception as e:
            print(f"警告：饱和温度计算失败 (P={pressure})")
            print(f"错误信息: {str(e)}")
            return 373.15  # 返回1个大气压下的近似值

    def calculate_single_phase_htc(self, i, j, phase):
        """
        计算单相对流换热系数
        
        参数：
        i, j: int, 网格位置
        phase: str, 'liquid' 或 'vapor'
        
        返回：
        float: 换热系数 (W/m²·K)
        """
        if phase == 'liquid':
            props = self.get_fluid_properties(self.pressure[i,j], self.T_l[i,j], 'liquid')
            velocity = self.liquid_velocity[i,j]
        else:
            props = self.get_fluid_properties(self.pressure[i,j], self.T_v[i,j], 'vapor')
            velocity = self.vapor_velocity[i,j]
        
        if props is None:
            return 0.0
        
        Re = props['density'] * velocity * self.channel_diameter / props['viscosity']
        Pr = props['viscosity'] * props.get('specific_heat', 4200) / props.get('thermal_conductivity', 0.6)
        
        # Dittus-Boelter correlation
        if Re > 10000 and 0.7 <= Pr <= 160:
            Nu = 0.023 * Re**0.8 * Pr**0.4
        else:
            Nu = 4.36  # 层流或过渡区域的简化处理
        
        h = Nu * props.get('thermal_conductivity', 0.6) / self.channel_diameter
        return h

    def check_convergence(self, old_void, old_vl, old_vv, old_p, old_t):
        """检查收敛性"""
        max_error = max(
            np.max(np.abs(self.void_fraction - old_void)),
            np.max(np.abs(self.liquid_velocity - old_vl)),
            np.max(np.abs(self.vapor_velocity - old_vv)),
            np.max(np.abs(self.pressure - old_p)) / self.inlet_pressure,
            np.max(np.abs(self.temperature - old_t)) / self.inlet_temp)
        
        return max_error

    def setup_staggered_grid(self):
        """设置交错网格"""
        # 标量网格（压力、温度、空泡率等）
        self.scalar_grid = np.zeros((self.n_axial, self.n_channels))
        
        # 轴向速度网格（错开半个网格）
        self.axial_velocity_grid = np.zeros((self.n_axial + 1, self.n_channels))
        
        # 横向速度网格（错开半个网格）
        self.lateral_velocity_grid = np.zeros((self.n_axial, self.n_channels + 1))
        
        # 网格中心点坐标
        self.z_centers = np.linspace(0, self.channel_length, self.n_axial)
        self.z_faces = np.linspace(-self.dz/2, self.channel_length + self.dz/2, self.n_axial + 1)
        
    def build_preconditioner(self):
        """构建基于物理过程的预处理器"""
        n = self.n_axial * self.n_channels
        M = lil_matrix((n, n))
        
        # 构建预处理器矩阵（仅考虑轴向相邻单元）
        for i in range(self.n_axial):
            for j in range(self.n_channels):
                idx = i * self.n_channels + j
                
                # 中心单元
                M[idx, idx] = 1.0
                
                # 轴向相邻单元
                if i > 0:
                    M[idx, (i-1) * self.n_channels + j] = -0.5
                if i < self.n_axial - 1:
                    M[idx, (i+1) * self.n_channels + j] = -0.5
                    
        return M.tocsr()
        
    def solve_pkm(self):
        """使用Picard Krylov方法求解方程组"""
        converged = False
        picard_iter = 0
        
        while not converged and picard_iter < self.max_picard_iterations:
            # 1. 计算热通量
            self.calculate_heat_flux()
            
            # 2. 求解动量方程
            self.solve_momentum_equations()
            
            # 3. 构建质量-能量方程组
            A, b = self.build_mass_energy_system()
            
            # 4. 构建预处理器
            M = self.build_preconditioner()
            
            # 5. 使用PGMRES求解线性方程组
            try:
                x, info = gmres(A, b, M=M, tol=self.krylov_tolerance, 
                              maxiter=self.max_krylov_iterations)
                
                if info > 0:
                    print(f"警告：GMRES在{info}次迭代后未收敛")
                    
                # 更新解
                self.update_solution(x)
                
                # 计算收敛性
                residual = self.calculate_residual()
                converged = residual < self.picard_tolerance
                
            except Exception as e:
                print(f"错误：GMRES求解失败 - {str(e)}")
                # 如果求解失败，减小时间步长重试
                self.time_step *= 0.5
                if self.time_step < 1e-6:
                    raise Exception("时间步长过小，计算终止")
                continue
                
            picard_iter += 1
            
        return converged, picard_iter
        
    def build_mass_energy_system(self):
        """构建质量-能量方程组的系数矩阵和右端向量"""
        n = self.n_axial * self.n_channels
        A = lil_matrix((n, n))
        b = np.zeros(n)
        
        for i in range(self.n_axial):
            for j in range(self.n_channels):
                idx = i * self.n_channels + j
                
                # 1. 质量守恒方程
                # 轴向质量流量
                if i > 0:
                    A[idx, idx-self.n_channels] = -self.liquid_velocity[i-1,j]
                A[idx, idx] = self.liquid_velocity[i,j]
                if i < self.n_axial - 1:
                    A[idx, idx+self.n_channels] = -self.liquid_velocity[i+1,j]
                
                # 横向质量流量
                for k in self.channel_connections[j+1]:
                    k_idx = i * self.n_channels + (k-1)
                    w_jk = self.calculate_cross_flow(i, j+1, k)
                    A[idx, k_idx] = w_jk
                    
                # 2. 能量守恒方程
                # 轴向能量输运
                if i > 0:
                    A[idx, idx-self.n_channels] += -self.liquid_velocity[i-1,j] * self.liquid_enthalpy[i-1,j]
                A[idx, idx] += self.liquid_velocity[i,j] * self.liquid_enthalpy[i,j]
                if i < self.n_axial - 1:
                    A[idx, idx+self.n_channels] += -self.liquid_velocity[i+1,j] * self.liquid_enthalpy[i+1,j]
                
                # 加热项
                q_wall = self.calculate_wall_heat_flux(i, j, 'liquid')
                b[idx] = q_wall
                
        return A.tocsr(), b
        
    def solve_momentum_equations(self):
        """求解动量方程"""
        # 轴向动量方程
        for i in range(1, self.n_axial):
            for j in range(self.n_channels):
                # 获取交错网格上的速度
                v_up = self.axial_velocity_grid[i-1, j]
                v_center = self.axial_velocity_grid[i, j]
                v_down = self.axial_velocity_grid[i+1, j] if i < self.n_axial-1 else v_center
                
                # 计算压力梯度
                dp_dz = (self.pressure[i,j] - self.pressure[i-1,j]) / self.dz
                
                # 计算摩擦力
                f_wall = self.calculate_wall_friction(i, j, 'liquid')
                
                # 更新轴向速度
                self.axial_velocity_grid[i,j] = v_center - self.time_step * (
                    v_center * (v_down - v_up)/(2*self.dz) +  # 对流项
                    dp_dz/self.get_fluid_properties(self.pressure[i,j], 
                        self.temperature[i,j], 'liquid')['density'] +  # 压力项
                    f_wall +  # 摩擦力
                    9.81  # 重力项
                )
        
        # 横向动量方程
        for i in range(self.n_axial):
            for j in range(self.n_channels-1):
                # 计算横向压力梯度
                dp_dy = self.calculate_lateral_pressure_gradient(i, j)
                
                # 计算横向摩擦力
                f_lat = self.calculate_wall_friction(i, j, 'liquid', 'lateral')
                
                # 更新横向速度
                self.lateral_velocity_grid[i,j] = self.lateral_velocity_grid[i,j] - self.time_step * (
                    dp_dy/self.get_fluid_properties(self.pressure[i,j], 
                        self.temperature[i,j], 'liquid')['density'] +
                    f_lat
                )
                
    def update_solution(self, x):
        """更新解向量"""
        for i in range(self.n_axial):
            for j in range(self.n_channels):
                idx = i * self.n_channels + j
                
                # 更新压力
                self.pressure[i,j] = x[idx]
                
                # 更新温度（通过压力和焓计算）
                h = x[idx + self.n_axial * self.n_channels]
                self.temperature[i,j] = self.get_temperature_from_enthalpy(
                    h, self.pressure[i,j], 'liquid')
                
    def calculate_heat_flux(self):
        """计算热通量"""
        # 计算热通量
        # ...

    def calculate_friction_factor(self, velocity, viscosity, flow_area, wetted_perimeter):
        """计算摩擦因子"""
        # 计算水力直径
        hydraulic_diameter = 4 * flow_area / wetted_perimeter
        
        # 计算雷诺数
        reynolds = abs(velocity) * hydraulic_diameter / viscosity
        
        # 层流区
        if reynolds < 2300:
            return 64 / reynolds
        # 湍流区
        else:
            return 0.316 * reynolds**(-0.25)

    def debug_output(self, time_step):
        """输出调试信息，用于验证16子通道计算"""
        try:
            print(f"\n====== 时间步 {time_step} 调试信息 ======")
            
            # 1. 检查质量守恒
            inlet_mass_flow = 0
            outlet_mass_flow = 0
            for j in range(self.n_channels):
                channel_type = self.channel_types[j+1]
                flow_area = self.channel_geometry[channel_type]['flow_area']
                
                # 入口质量流量
                rho_in = self.get_fluid_properties(self.pressure[0,j], self.T_l[0,j], 'liquid')['density']
                inlet_mass_flow += rho_in * self.liquid_velocity[0,j] * flow_area
                
                # 出口质量流量
                rho_l_out = self.get_fluid_properties(self.pressure[-1,j], self.T_l[-1,j], 'liquid')['density']
                rho_v_out = self.get_fluid_properties(self.pressure[-1,j], self.T_v[-1,j], 'vapor')['density']
                outlet_mass_flow += (rho_l_out * (1-self.void_fraction[-1,j]) * self.liquid_velocity[-1,j] + 
                                    rho_v_out * self.void_fraction[-1,j] * self.vapor_velocity[-1,j]) * flow_area
            
            mass_error = abs(inlet_mass_flow - outlet_mass_flow) / inlet_mass_flow
            print(f"质量守恒误差: {mass_error:.6f}")
            
            # 2. 检查能量守恒
            total_power_in = 0
            total_power_out = 0
            total_heat_added = 0
            
            for j in range(self.n_channels):
                channel_type = self.channel_types[j+1]
                flow_area = self.channel_geometry[channel_type]['flow_area']
                heated_perimeter = self.channel_geometry[channel_type]['heated_perimeter']
                
                # 入口能量
                h_in = self.h_l[0,j]
                total_power_in += inlet_mass_flow * h_in / self.n_channels
                
                # 出口能量
                h_l_out = self.h_l[-1,j]
                h_v_out = self.h_v[-1,j]
                total_power_out += outlet_mass_flow * (
                    (1-self.void_fraction[-1,j]) * h_l_out + 
                    self.void_fraction[-1,j] * h_v_out
                ) / self.n_channels
                
                # 加热功率
                for i in range(self.n_axial):
                    total_heat_added += self.q_wl[i,j] * self.dz
            
            energy_error = abs(total_power_out - total_power_in - total_heat_added) / total_heat_added
            print(f"能量守恒误差: {energy_error:.6f}")
            
            # 3. 检查各子通道类型的平均值
            for type_name in ['corner', 'edge', 'center']:
                void_sum = 0
                temp_sum = 0
                count = 0
                
                for j in range(self.n_channels):
                    if self.channel_types[j+1] == type_name:
                        void_sum += np.mean(self.void_fraction[:,j])
                        temp_sum += np.mean(self.T_l[:,j])
                        count += 1
                
                if count > 0:
                    print(f"\n{type_name}型子通道平均值:")
                    print(f"平均空泡率: {void_sum/count:.4f}")
                    print(f"平均温度: {temp_sum/count:.2f} K")
            
            # 4. 检查数值稳定性
            max_courant = 0
            for j in range(self.n_channels):
                channel_type = self.channel_types[j+1]
                flow_area = self.channel_geometry[channel_type]['flow_area']
                
                # 计算Courant数
                max_velocity = max(np.max(abs(self.liquid_velocity[:,j])), 
                                 np.max(abs(self.vapor_velocity[:,j])))
                courant = max_velocity * self.dt / self.dz
                max_courant = max(max_courant, courant)
            
            print(f"\n最Courant数: {max_courant:.4f}")
            if max_courant > 1:
                print("警告：Courant数大于1，可能存在数值不稳定")
                
        except Exception as e:
            print(f"尝试输出时出错: {str(e)}")

    def verify_channel_connections(self):
        """验证子通道连接的���确性"""
        try:
            print("\n====== 验证子通道连接 ======")
            
            # 检查连接的对称性
            for j in range(1, self.n_channels + 1):
                for k in self.channel_connections[j]:
                    if j not in self.channel_connections[k]:
                        print(f"错误：子通道{j}和{k}的连接不对称")
                        
            # 检查连接的物理合理性
            for j in range(1, self.n_channels + 1):
                connections = self.channel_connections[j]
                channel_type = self.channel_types[j]
                
                # 检查连接数量
                if channel_type == 'corner' and len(connections) != 2:
                    print(f"警告：角通道{j}的连接数量不正确")
                elif channel_type == 'edge' and len(connections) != 3:
                    print(f"警告：边通道{j}的连接数量不正确")
                elif channel_type == 'center' and len(connections) != 4:
                    print(f"警告：中心通道{j}的连接数量不正确")
                    
            print("子通道连接验证完成")
            
        except Exception as e:
            print(f"验证子通道连接时出错: {str(e)}")

    def check_numerical_stability(self):
        """检查数值稳定性并调整时间步长"""
        try:
            # 1. 检查物理量的变化率
            max_void_change = np.max(np.abs(np.diff(self.void_fraction, axis=0)))
            max_velocity_change = max(
                np.max(np.abs(np.diff(self.liquid_velocity, axis=0))),
                np.max(np.abs(np.diff(self.vapor_velocity, axis=0)))
            )
            max_pressure_change = np.max(np.abs(np.diff(self.pressure, axis=0))) / np.mean(self.pressure)
            max_temp_change = max(
                np.max(np.abs(np.diff(self.T_l, axis=0))) / np.mean(self.T_l),
                np.max(np.abs(np.diff(self.T_v, axis=0))) / np.mean(self.T_v))
            
            
            # 2. 设置阈值
            thresholds = {
                'void': 0.1,
                'velocity': 0.2,
                'pressure': 0.1,
                'temperature': 0.1
            }
            
            # 3. 根据变化率调整时间步长
            dt_factor = 1.0
            
            if max_void_change > thresholds['void']:
                dt_factor = min(dt_factor, 0.5)
                logging.warning(f"空泡率变化过大: {max_void_change:.4f}")
                
            if max_velocity_change > thresholds['velocity']:
                dt_factor = min(dt_factor, 0.5)
                logging.warning(f"速度变化过大: {max_velocity_change:.4f}")
                
            if max_pressure_change > thresholds['pressure']:
                dt_factor = min(dt_factor, 0.5)
                logging.warning(f"压力变化过大: {max_pressure_change:.4f}")
                
            if max_temp_change > thresholds['temperature']:
                dt_factor = min(dt_factor, 0.5)
                logging.warning(f"温度变化过大: {max_temp_change:.4f}")
                
            # 4. 检查物理约束
            if (np.any(self.void_fraction < 0) or np.any(self.void_fraction > 1) or
                np.any(self.pressure < 0) or
                np.any(self.T_l < 273.15) or np.any(self.T_v < 273.15)):
                dt_factor = 0.1
                logging.error("物理量超出合理范围")
                
            # 5. 限制时间步长的变化范围
            self.dt = np.clip(self.dt * dt_factor, 1e-6, 1e-2)
            
            if dt_factor < 1.0:
                logging.warning(f"时间步长已调整为: {self.dt:.2e}")
                
            return dt_factor < 1.0
            
        except Exception as e:
            logging.error(f"检查数值稳定性时出错: {str(e)}")
            return True  # 返回True表示需要调整时间步长

    def calculate_liquid_htc(self, i, j):
        """
        计算液相换热系数
        
        使用Dittus-Boelter关联式计算单相强制对流换热系数。
        考虑不同类型子通道的几何特征。
        
        参数：
            i (int): 轴向位置索引
            j (int): 子通道编号
            
        返回：
            float: 液相换热系数 (W/m²·K)
            
        注意：
            - 使用Dittus-Boelter关联式: Nu = 0.023 * Re^0.8 * Pr^0.4
            - 考虑不同类型子通道的水力直径
            - 包含异常处理机制
        """
        try:
            # 获取子通道类型和几何参数
            channel_type = self.channel_types[j]
            flow_area = self.channel_geometry[channel_type]['flow_area']
            wetted_perimeter = self.channel_geometry[channel_type]['wetted_perimeter']
            
            # 计算水力直径
            d_h = 4.0 * flow_area / wetted_perimeter
            
            # 获取流体物性
            p = self.pressure[i, j-1]
            t = self.temperature[i, j-1]
            props = self.get_fluid_properties(p, t, 'liquid')
            
            if props is None:
                return 0.0
                
            rho = props['density']
            mu = props['viscosity']
            
            # 计算雷诺数
            v = self.liquid_velocity[i, j-1]
            re = rho * abs(v) * d_h / mu
            
            # 计算普朗特数 (假设为常数)
            pr = 1.0
            
            # 使用Dittus-Boelter关联式
            if re > 2300:  # 湍流
                nu = 0.023 * re**0.8 * pr**0.4
            else:  # 层流
                nu = 3.66
                
            # 计算导热系数 (假设为常数)
            k = 0.6  # W/m·K
            
            # 计算换热系数
            h = nu * k / d_h
            
            return h
            
        except Exception as e:
            print(f"计算液相换热系数时出错: {str(e)}")
            return 0.0
            
    def calculate_vapor_htc(self, i, j):
        """
        计算气相换热系数
        
        使用改进的Dittus-Boelter关联式计算气相的换热系数。
        考虑气相的特殊流动特性。
        
        参数：
            i (int): 轴向位置索引
            j (int): 子通道编号
            
        返回：
            float: 气相换热系数 (W/m²·K)
            
        注意：
            - 使用修正的关联式: Nu = 0.018 * Re^0.8 * Pr^0.4
            - 系数0.018考虑了气相流动的特殊性
            - 包含对低雷诺数情况的处理
        """
        try:
            # 获取子通道几何参数
            channel_type = self.channel_types[j]
            flow_area = self.channel_geometry[channel_type]['flow_area']
            wetted_perimeter = self.channel_geometry[channel_type]['wetted_perimeter']
            
            # 计算水力直径
            d_h = 4.0 * flow_area / wetted_perimeter
            
            # 获取流体物性
            p = self.pressure[i, j-1]
            t = self.temperature[i, j-1]
            props = self.get_fluid_properties(p, t, 'vapor')
            
            if props is None:
                return 0.0
                
            rho = props['density']
            mu = props['viscosity']
            
            # 计算雷诺数
            v = self.vapor_velocity[i, j-1]
            re = rho * abs(v) * d_h / mu
            
            # 计算普朗特数 (假设为常数)
            pr = 1.0
            
            # 使用修正的关联式
            if re > 2300:  # 湍流
                nu = 0.018 * re**0.8 * pr**0.4  # 系数改为0.018
            else:  # 层流
                nu = 3.66
                
            # 计算导热系数 (假设为常数)
            k = 0.025  # W/m·K (气相导热系数较小)
            
            # 计算换热系数
            h = nu * k / d_h
            
            return h
            
        except Exception as e:
            print(f"计算气相换热系数时出错: {str(e)}")
            return 0.0
            
    def calculate_wall_temperature(self, i, j):
        """
        计算壁面温度
        
        基于热流密度和换热系数计算壁面温度。
        考虑两相流动的换热特性。
        
        参数：
            i (int): 轴向位置索引
            j (int): 子通道编号
            
        返回：
            float: 壁面温度 (K)
            
        注意：
            - 考虑液相和气相的共同换热效果
            - 使用加权平均的方法处理两相换热
            - 包含对极限情况的处理
        """
        try:
            # 获取流体温度
            t_f = self.temperature[i, j-1]
            
            # 获取热流密度
            q_wl = self.q_wl[i, j-1]  # 液相热流密度
            q_wv = self.q_wv[i, j-1]  # 气相热流密度
            
            # 计算换热系数
            h_l = self.calculate_liquid_htc(i, j)
            h_v = self.calculate_vapor_htc(i, j)
            
            # 获取空泡率
            alpha = self.void_fraction[i, j-1]
            
            # 计算总的换热系数
            h_total = (1 - alpha) * h_l + alpha * h_v
            
            # 计算总的热流密度
            q_total = q_wl + q_wv
            
            # 计算壁面温度
            if h_total > 0:
                t_w = t_f + q_total / h_total
            else:
                t_w = t_f + 10  # 默认温差10K
                
            return t_w
            
        except Exception as e:
            print(f"计算壁面温度时出错: {str(e)}")
            return self.temperature[i, j-1] + 10  # 返回一个安全的估计值

    def calculate_condensation_heat_flux(self, T_w, T_sat, P):
        """计算冷凝热流密度"""
        try:
            # 获取气相物性
            props_v = self.get_fluid_properties(P, T_sat, 'vapor')
            if props_v is None or 'thermal_conductivity' not in props_v:
                print(f"警告：流体物性数据不完整 (P={P}, T_sat={T_sat})")
                return 0.0
            
            # 计算冷凝换热系数
            h_cond = self.calculate_single_phase_htc(T_w, P, 'vapor')  # 假设使用气相的单相换热系数
            
            # 计��冷凝热流密度
            q_cond = h_cond * (T_w - T_sat)
            
            return max(0, q_cond)  # 确保热流密度非负
            
        except Exception as e:
            print(f"计算冷凝热流密度时出错: {str(e)}")
            return 0.0

    def calculate_boiling_heat_flux(self, T_w, T_sat, P):
        """计算汽化热流密度"""
        try:
            # 获取液相物性
            props_l = self.get_fluid_properties(P, T_sat, 'liquid')
            if props_l is None or 'thermal_conductivity' not in props_l:
                print(f"警告：流体物性数据不完整 (P={P}, T_sat={T_sat})")
                return 0.0
            
            # 计算汽化换热系数
            h_boil = self.calculate_single_phase_htc(T_w, P, 'liquid')  # 假设使用液相的单相换热系数
            
            # 计算汽化热流密度
            q_boiling = h_boil * (T_w - T_sat)
            
            return max(0, q_boiling)  # 确保热流密度非负
            
        except Exception as e:
            print(f"计算汽化热流密度时出错: {str(e)}")
            return 0.0

    def get_temperature_from_enthalpy(self, enthalpy, pressure, phase):
        """
        根据给定的焓和压力计算温度
        
        参数：
            enthalpy (float): 比焓 (J/kg)
            pressure (float): 压力 (Pa)
            phase (str): 相态 ('liquid' 或 'vapor')
        
        返回：
            float: 温度 (K)
        """
        try:
            from CoolProp.CoolProp import PropsSI

            # 检查焓和压力的有效性
            if enthalpy <= 0 or pressure <= 0:
                raise ValueError("焓和压力必须为正值")

            # 使用PropsSI检查给定压力下的饱和焓
            T_sat = self.get_saturation_temperature(pressure)
            h_sat = PropsSI('H', 'P', pressure, 'Q', 0, 'Water')  # 饱和液体焓

            if enthalpy < h_sat:
                print(f"警告：给定焓 {enthalpy} 小于饱和液体焓 {h_sat}，可能不合理。")
                return T_sat  # 返回饱和温度

            def objective(T):
                if phase == 'liquid':
                    return PropsSI('H', 'P', pressure, 'T', T, 'Water') - enthalpy
                elif phase == 'vapor':
                    return PropsSI('H', 'P', pressure, 'T', T, 'Water') - enthalpy
                else:
                    raise ValueError("无效的相态参数，需为 'liquid' 或 'vapor'")

            # 初始猜测温度
            T_initial = self.get_saturation_temperature(pressure)
            T_solution = opt.newton(objective, T_initial)
            return T_solution

        except Exception as e:
            print(f"警告：根据焓和压力计算温度失败 (H={enthalpy}, P={pressure})")
            print(f"错误信息: {str(e)}")
            return 373.15  # 返回默认温度

    def calculate_interface_tension(self, i, j):
        """计算界面张力效应"""
        try:
            # 获取局部条件
            alpha = self.void_fraction[i,j]
            d_bubble = 0.001  # 假设的气泡直径 (m)
            
            # 计算界面曲率
            if i > 0 and i < self.n_axial-1:
                d2alpha_dz2 = (self.void_fraction[i+1,j] - 2*self.void_fraction[i,j] + 
                             self.void_fraction[i-1,j]) / (self.dz**2)
            else:
                d2alpha_dz2 = 0
                
            # 计算界面张力
            F_st = -self.surface_tension * d2alpha_dz2
            
            # 考虑壁面接触角的影响
            theta = np.radians(self.contact_angle)
            F_st *= np.cos(theta)
            
            return F_st
            
        except Exception as e:
            print(f"计算界面张力时出错: {str(e)}")
            return 0.0

    def solve_k_epsilon(self, i, j):
        """求解k-ε湍流模型方程"""
        try:
            # 获取局部流动参数
            u = self.liquid_velocity[i,j]
            rho = self.get_fluid_properties(self.pressure[i,j], 
                                          self.temperature[i,j], 'liquid')['density']
            mu = self.get_fluid_properties(self.pressure[i,j], 
                                         self.temperature[i,j], 'liquid')['viscosity']
            
            # 计算速度梯度
            if i > 0 and i < self.n_axial-1:
                du_dy = (self.liquid_velocity[i,j+1] - self.liquid_velocity[i,j-1]) / (2*self.dz)
                
                # 计算湍流产生项
                P_k = mu * (du_dy**2)
                
                # 求解k方程
                self.k[i,j] += self.dt * (
                    P_k - self.epsilon[i,j] +
                    self.calculate_k_diffusion(i, j)
                )
                
                # 求解ε方程
                self.epsilon[i,j] += self.dt * (
                    self.C_1 * self.epsilon[i,j] * P_k / self.k[i,j] -
                    self.C_2 * self.epsilon[i,j]**2 / self.k[i,j] +
                    self.calculate_epsilon_diffusion(i, j)
                )
                
                # 计算湍流粘度
                mu_t = self.C_mu * rho * self.k[i,j]**2 / self.epsilon[i,j]
                
                return mu_t
                
        except Exception as e:
            print(f"求解k-ε模型时出错: {str(e)}")
            return 0.0

    def solve_heat_coupling(self, i, j):
        """求解热耦合效应"""
        try:
            # 计算燃料到包壳的导热
            q_fuel_clad = self.k_fuel * (self.T_fuel[i,j] - self.T_clad[i,j]) / self.gap_conductance
            
            # 计算包壳到冷却剂的导热
            q_clad_cool = self.k_clad * (self.T_clad[i,j] - self.temperature[i,j])
            
            # 更新温度场
            self.T_fuel[i,j] += self.dt * (
                self.calculate_fuel_heat_generation(i, j) - q_fuel_clad
            ) / (self.rho_fuel * self.cp_fuel)
            
            self.T_clad[i,j] += self.dt * (
                q_fuel_clad - q_clad_cool
            ) / (self.rho_clad * self.cp_clad)
            
            return q_clad_cool
            
        except Exception as e:  # 添加异常处理
            print(f"求解热耦合效应时出错: {str(e)}")
            return 0.0  # 返一个默认值以防出错

    def update_values(self, prev_void_fraction, prev_liquid_velocity, prev_vapor_velocity, time_step):
        """更新值"""
        new_void = np.zeros_like(prev_void_fraction)
        new_vl = np.zeros_like(prev_liquid_velocity)
        new_vv = np.zeros_like(prev_vapor_velocity)

        for i in range(new_void.shape[0]):
            for j in range(new_void.shape[1]):
                d_alpha_dx = 0
                d_vl_dx = 0
                d_vv_dx = 0
                dp_dx = 0

                # 更新值
                new_void[i,j] = prev_void_fraction[i,j] - time_step * (
                    vapor_velocity[i,j] * d_alpha_dx
                )

                new_vl[i,j] = prev_liquid_velocity[i,j] - time_step * (
                    liquid_velocity[i,j] * d_vl_dx +
                    dp_dx +
                    9.81
                )

                new_vv[i,j] = prev_vapor_velocity[i,j] - time_step * (
                    vapor_velocity[i,j] * d_vv_dx +
                    dp_dx +
                    9.81
                )

                # 限制物理范围
                new_void[i,j] = min(max(new_void[i,j], 0.0), 0.9999)
                new_vl[i,j] = max(new_vl[i,j], 0.0)
                new_vv[i,j] = max(new_vv[i,j], 0.0)

        return new_void, new_vl, new_vv
            
        

@njit(parallel=True)
def compute_energy_points_numba(n_axial, n_channels, temperature, prev_temperature,
                              liquid_velocity, pressure, time_step, dz):
    """使用Numba优化的能量点计算"""
    new_temp = np.zeros_like(temperature)
    
    for i in prange(1, n_axial):
        for j in prange(n_channels):
            # 计算空间导数项
            if i > 0:
                d_T_dx = (temperature[i,j] - temperature[i-1,j]) / dz
            else:
                d_T_dx = 0
            
            # 更新温度
            new_temp[i,j] = prev_temperature[i,j] - time_step * (
                liquid_velocity[i,j] * d_T_dx
            )
            
    return new_temp

@njit
def calculate_wall_friction_numba(velocity, diameter, density, viscosity):
    """使用Numba优化的壁面摩擦力计算"""
    # 计算Reynolds数
    Re = density * abs(velocity) * diameter / viscosity
    
    # 计算摩擦系数
    if Re < 2300:  # 层流
        f = 64 / Re
    else:  # 湍流
        f = 0.316 * Re**(-0.25)
        
    # 计算壁面摩擦力
    tau_w = 0.5 * f * density * velocity * abs(velocity)
    
    return tau_w

@njit(parallel=True)
def calculate_residual_numba(n_axial, n_channels, liquid_velocity, liquid_enthalpy, dz):
    """使用Numba优化的残差计算"""
    mass_residual = np.zeros((n_axial, n_channels))
    energy_residual = np.zeros((n_axial, n_channels))
    
    for i in prange(1, n_axial):
        for j in prange(n_channels):
            # 轴向质量流量差
            dm_dz = (liquid_velocity[i,j] - liquid_velocity[i-1,j]) / dz
            mass_residual[i,j] = dm_dz
            
            # 轴向能量输运
            dh_dz = (liquid_enthalpy[i,j] - liquid_enthalpy[i-1,j]) / dz
            energy_residual[i,j] = dh_dz
            
    return max(np.max(np.abs(mass_residual)), np.max(np.abs(energy_residual)))

@njit(parallel=True)
def solve_two_fluid_equations_numba(n_axial, n_channels, 
                                  alpha_l, alpha_v, rho_l, rho_v,
                                  u_l, u_v, v_l, v_v, w_l, w_v,
                                  P, T_l, T_v, T_c, h_l, h_v,
                                  dx, dy, dz, dt):
    """使用Numba优化的双流体模型求解器"""
    # 计算CFL时间步长
    dt_cfl = calculate_cfl_timestep(u_l, u_v, v_l, v_v, dx, dy, dz)
    
    # 使用较小的时间步长
    dt = min(dt, dt_cfl)
    
    # 初始化新的场值
    new_alpha_l = np.zeros_like(alpha_l)
    new_alpha_v = np.zeros_like(alpha_v)
    new_u_l = np.zeros_like(u_l)
    new_u_v = np.zeros_like(u_v)
    new_v_l = np.zeros_like(v_l)
    new_v_v = np.zeros_like(v_v)
    new_h_l = np.zeros_like(h_l)
    new_h_v = np.zeros_like(h_v)
    
    for i in prange(1, n_axial-1):
        for j in prange(1, n_channels-1):
            # ====== 1. 计算相变质量流量 ======
            T_sat = calculate_saturation_temperature_numba(P[i,j])
            
            # 过冷沸腾
            if T_l[i,j] < T_sat and T_c[i,j] > T_sat:
                h_sc = calculate_subcooled_boiling_htc_numba(T_c[i,j], T_sat, T_l[i,j], P[i,j])
                Q_sc = h_sc * (T_c[i,j] - T_l[i,j])
                gamma_sc = Q_sc / calculate_latent_heat_numba(P[i,j])
            else:
                gamma_sc = 0.0
                
            # 饱和核态沸腾
            if T_l[i,j] >= T_sat and T_c[i,j] > T_sat:
                h_sh = calculate_nucleate_boiling_htc_numba(T_c[i,j], T_sat, P[i,j])
                Q_sh = h_sh * (T_c[i,j] - T_sat)
                gamma_sh = Q_sh / calculate_latent_heat_numba(P[i,j])
            else:
                gamma_sh = 0.0
                
            # 净相变质量流量
            gamma = gamma_sh - gamma_sc
            
            # ====== 2. 计算湍流混合项 ======
            # 液相湍流混合
            T_ml = calculate_turbulent_mixing_numba(alpha_l[i,j], alpha_l[i,j+1], 
                                                  rho_l[i,j], dy)
            # 气相湍流混合
            T_mv = calculate_turbulent_mixing_numba(alpha_v[i,j], alpha_v[i,j+1], 
                                                  rho_v[i,j], dy)
                    
            # ====== 3. 计算界面力 ======
            F_lv = calculate_interfacial_force_numba(alpha_l[i,j], rho_l[i,j],
                                                   u_l[i,j], u_v[i,j],
                                                   v_l[i,j], v_v[i,j])
            
            # ====== 4. 液相质量守恒方程 ======
            d_alpha_l_dt = -(1/rho_l[i,j]) * (
                # 轴向质量通量导数
                (alpha_l[i+1,j]*rho_l[i+1,j]*u_l[i+1,j] - 
                 alpha_l[i-1,j]*rho_l[i-1,j]*u_l[i-1,j])/(2*dz) +
                # 横向质量通量导数
                (alpha_l[i,j+1]*rho_l[i,j+1]*v_l[i,j+1] - 
                 alpha_l[i,j-1]*rho_l[i,j-1]*v_l[i,j-1])/(2*dy)
            ) - gamma/rho_l[i,j] + T_ml
            
            # ====== 5. 气相质量守恒方程 ======
            d_alpha_v_dt = -(1/rho_v[i,j]) * (
                # 轴向量通量导数
                (alpha_v[i+1,j]*rho_v[i+1,j]*u_v[i+1,j] - 
                 alpha_v[i-1,j]*rho_v[i-1,j]*u_v[i-1,j])/(2*dz) +
                # 横向质量通量导数
                (alpha_v[i,j+1]*rho_v[i,j+1]*v_v[i,j+1] - 
                 alpha_v[i,j-1]*rho_v[i,j-1]*v_v[i,j-1])/(2*dy)
            ) + gamma/rho_v[i,j] + T_mv
            
            # ====== 6. 液相轴向动量方程 ======
            # 壁面摩擦力
            F_wl = calculate_wall_friction_numba(u_l[i,j], 0.01, rho_l[i,j], 1e-3)
            
            d_ul_dt = -(
                # 对流项
                u_l[i,j] * (u_l[i+1,j] - u_l[i-1,j])/(2*dz) +
                v_l[i,j] * (u_l[i,j+1] - u_l[i,j-1])/(2*dy) +
                # 压力梯度
                (1/rho_l[i,j]) * (P[i+1,j] - P[i-1,j])/(2*dz) +
                # 重力项
                9.81 +
                # 壁面摩擦力
                F_wl/rho_l[i,j] -
                # 界面力
                F_lv/rho_l[i,j] -
                # 相变动量
                gamma * u_l[i,j]/rho_l[i,j]
            )
            
            # ====== 7. 气相轴向动量方程 ======
            # 壁面摩擦力
            F_wv = calculate_wall_friction_numba(u_v[i,j], 0.01, rho_v[i,j], 1e-5)
            
            d_uv_dt = -(
                # 对流项
                u_v[i,j] * (u_v[i+1,j] - u_v[i-1,j])/(2*dz) +
                v_v[i,j] * (u_v[i,j+1] - u_v[i,j-1])/(2*dy) +
                # 压力梯度
                (1/rho_v[i,j]) * (P[i+1,j] - P[i-1,j])/(2*dz) +
                # 重力项
                9.81 +
                # 壁面摩擦力
                F_wv/rho_v[i,j] +
                # 界面力
                F_lv/rho_v[i,j] +
                # 相变动量
                gamma * u_v[i,j]/rho_v[i,j]
            )
            
            # ====== 8. 液相横向动量方程 ======
            d_vl_dt = -(
                # 对流项
                u_l[i,j] * (v_l[i+1,j] - v_l[i-1,j])/(2*dz) +
                v_l[i,j] * (v_l[i,j+1] - v_l[i,j-1])/(2*dy) -
                # 界面力
                F_lv/rho_l[i,j] -
                # 相变动量
                gamma * v_l[i,j]/rho_l[i,j]
            )
            
            # ====== 9. 气相横向动量方程 ======
            d_vv_dt = -(
                # 对流
                u_v[i,j] * (v_v[i+1,j] - v_v[i-1,j])/(2*dz) +
                v_v[i,j] * (v_v[i,j+1] - v_v[i,j-1])/(2*dy) +
                # 压力梯度
                (1/rho_v[i,j]) * (P[i,j+1] - P[i,j-1])/(2*dy) +
                # 界面力
                F_lv/rho_v[i,j] +
                # 相变动量
                gamma * v_v[i,j]/rho_v[i,j]
            )
            
            # ====== 10. 液相能量方程 ======
            # 壁面换热
            q_wl = calculate_wall_heat_flux_numba(T_c[i,j], T_l[i,j], h_l[i,j])
            
            # 获取饱和焓
            try:
                h_sat_l = PropsSI('H', 'P', P[i,j], 'Q', 0, 'Water')
                h_sat_v = PropsSI('H', 'P', P[i,j], 'Q', 1, 'Water')
            except Exception as e:
                logging.error(f"计算饱和焓失败: {str(e)}")
                continue
            
            # 计算液相焓变化
            d_hl_dt = -(
                u_l[i,j] * (h_l[i+1,j] - h_l[i-1,j])/(2*dz) +  # 对流
                q_wl/(rho_l[i,j] * (1 - alpha_l[i,j])) +  # 壁面换热
                gamma * (h_sat_v - h_l[i,j])/(rho_l[i,j] * (1 - alpha_l[i,j]))  # 相变
            )

            # 计算气相焓变化
            d_hv_dt = -(
                u_v[i,j] * (h_v[i+1,j] - h_v[i-1,j])/(2*dz) +  # 对流
                q_wv/(rho_v[i,j] * alpha_v[i,j]) +  # 壁面换热
                gamma * (h_sat_l - h_v[i,j])/(rho_v[i,j] * alpha_v[i,j])  # 相变
            )

            # 更新焓值并限制在物理范围内
            h_l[i,j] = np.clip(h_l[i,j] + dt * d_hl_dt, h_sat_l * 0.8, h_sat_v)
            h_v[i,j] = np.clip(h_v[i,j] + dt * d_hv_dt, h_sat_l, h_sat_v * 1.2)

            # 更新温度
            try:
                T_l[i,j] = PropsSI('T', 'P', P[i,j], 'H', h_l[i,j], 'Water')
                T_v[i,j] = PropsSI('T', 'P', P[i,j], 'H', h_v[i,j], 'Water')
            except Exception as e:
                logging.warning(f"温度计算失败: {str(e)}")
                T_l[i,j] = T_l[i-1,j]
                T_v[i,j] = T_v[i-1,j]

            # 限制温度在物理范围内
            T_l[i,j] = np.clip(T_l[i,j], 273.15, 647.0)
            T_v[i,j] = np.clip(T_v[i,j], 273.15, 647.0)
            
            # ====== 11. 气相能量方程 ======
            # 壁面换热
            q_wv = calculate_wall_heat_flux_numba(T_c[i,j], T_v[i,j], h_v[i,j])
            
            d_hv_dt = -(
                # 对流项
                u_v[i,j] * (h_v[i+1,j] - h_v[i-1,j])/(2*dz) +
                v_v[i,j] * (h_v[i,j+1] - h_v[i,j-1])/(2*dy) +
                # 压力功
                alpha_v[i,j] * (P[i,j] - P[i-1,j])/dt +
                # 壁面换热
                q_wv/(rho_v[i,j] * alpha_v[i,j]) +
                # 相变热量
                gamma * h_v[i,j]/(rho_v[i,j] * alpha_v[i,j]) +
                # 添加燃料热生成和径向传导（气相部分）
                (q_fuel * alpha_v[i,j] + q_rad)/(rho_v[i,j] * alpha_v[i,j])
            )
            
            # ====== 12. 时间推进 ======
            new_alpha_l[i,j] = alpha_l[i,j] + dt * d_alpha_l_dt
            new_alpha_v[i,j] = alpha_v[i,j] + dt * d_alpha_v_dt
            new_u_l[i,j] = u_l[i,j] + dt * d_ul_dt
            new_u_v[i,j] = u_v[i,j] + dt * d_uv_dt
            new_v_l[i,j] = v_l[i,j] + dt * d_vl_dt
            new_v_v[i,j] = v_v[i,j] + dt * d_vv_dt
            new_h_l[i,j] = h_l[i,j] + dt * d_hl_dt
            new_h_v[i,j] = h_v[i,j] + dt * d_hv_dt
            
            # 限制物理范围
            new_alpha_l[i,j] = min(max(new_alpha_l[i,j], 0.0), 1.0)
            new_alpha_v[i,j] = 1.0 - new_alpha_l[i,j]
            
    # 在返回结果之前检查数值稳定性
    dt_new, needs_restart = check_numerical_stability(
        new_alpha_l, new_alpha_v, new_u_l, new_u_v, new_v_l, new_v_v,
        P, T_l, T_v, rho_l, rho_v, h_l, h_v, dt
    )
    
    if needs_restart:
        # 如果需要重新计算，返回特殊标记
        return (None, None, None, None, None, None, None, dt_new, True)
    
    return (new_alpha_l, new_alpha_v, new_u_l, new_u_v, 
            new_v_l, new_v_v, new_h_l, new_h_v, dt_new, False)

@njit
def calculate_turbulent_mixing_numba(alpha_i, alpha_j, rho, dy):
    """计算湍流混合项"""
    C_e = 0.01  # 涡流扩散率
    L_tur = dy  # 湍流混合长度
    alpha_eq_i = 0.5  # 平衡空泡分数
    
    return C_e * dy * L_tur * rho * (alpha_j - alpha_i - (alpha_j - alpha_eq_i))

@njit
def calculate_interfacial_force_numba(alpha_l, rho_l, u_l, u_v, v_l, v_v):
    """计算界面力"""
    C_i = 0.44  # 界面阻力系数
    d_b = 0.001  # 气泡直径
    
    # 计算相对速度
    u_r = u_v - u_l
    v_r = v_v - v_l
    V_r = np.sqrt(u_r**2 + v_r**2)
    
    # 计算界面面积密度
    a_i = 6 * (1-alpha_l) / d_b
    
    return 0.5 * C_i * rho_l * V_r * V_r * a_i

@njit
def calculate_wall_heat_flux_numba(T_w, T_f, h):
    """计算壁面热流密度"""
    return h * (T_w - T_f)

@njit
def calculate_saturation_temperature_numba(P):
    """计算饱和温度"""
    return 373.15 + 0.0001 * (P - 1.013e5)  # 简单的饱和温度计算

@njit
def calculate_latent_heat_numba(P):
    """计算汽化潜热"""
    return 2.257e6  # 常压下水汽化潜热

@njit
def calculate_cfl_timestep(u_l, u_v, v_l, v_v, dx, dy, dz):
    """计算基于CFL条件的时步长"""
    # 计算最大速度
    max_velocity = max(
        np.max(np.abs(u_l)),
        np.max(np.abs(u_v)),
        np.max(np.abs(v_l)),
        np.max(np.abs(v_v))
        )
    
    
    # CFL条件，取最小网格间距
    min_dx = min(dx, dy, dz)
    
    # CFL数设为0.5以保证稳定性
    cfl = 0.5
    
    # 计算时间步长
    dt = cfl * min_dx / (max_velocity + 1e-10)
    
    return dt

@njit
def check_numerical_stability(alpha_l, alpha_v, u_l, u_v, v_l, v_v, P, T_l, T_v,
                            rho_l, rho_v, h_l, h_v, dt_old):
    """检查数值稳定性并调整时间步长"""
    # 1. 检查物理量的变化率
    max_alpha_change = max(
        np.max(np.abs(np.diff(alpha_l))),
        np.max(np.abs(np.diff(alpha_v)))
    )
    
    max_velocity_change = max(
        np.max(np.abs(np.diff(u_l))),
        np.max(np.abs(np.diff(u_v))),
        np.max(np.abs(np.diff(v_l))),
        np.max(np.abs(np.diff(v_v))))
    
    
    max_pressure_change = np.max(np.abs(np.diff(P))) / np.mean(P)
    
    max_temp_change = max(
        np.max(np.abs(np.diff(T_l))) / np.mean(T_l),
        np.max(np.abs(np.diff(T_v))) / np.mean(T_v)
        )
    
    
    # 2. 设置阈值
    alpha_threshold = 0.1
    velocity_threshold = 0.2
    pressure_threshold = 0.1
    temp_threshold = 0.1
    
    # 3. 根据变化率调整时间步长
    dt_factor = 1.0
    
    if max_alpha_change > alpha_threshold:
        dt_factor = min(dt_factor, 0.5)
    if max_velocity_change > velocity_threshold:
        dt_factor = min(dt_factor, 0.5)
    if max_pressure_change > pressure_threshold:
        dt_factor = min(dt_factor, 0.5)
    if max_temp_change > temp_threshold:
        dt_factor = min(dt_factor, 0.5)
        
    # 如果所有变化率都很小，可以适当增加时间步长
    if (max_alpha_change < alpha_threshold/2 and
        max_velocity_change < velocity_threshold/2 and
        max_pressure_change < pressure_threshold/2 and
        max_temp_change < temp_threshold/2):
        dt_factor = min(1.2, 1.0/dt_factor)
        
    # 4. 检查物理约束
    if (np.any(alpha_l < 0) or np.any(alpha_l > 1) or
        np.any(P < 0) or
        np.any(T_l < 273) or np.any(T_v < 273) or
        np.any(rho_l < 0) or np.any(rho_v < 0)):
        dt_factor = 0.1  # 严重违反物理约束，大幅减小时间步长
        
    # 5. 限制时间步长的变化范围
    dt_new = dt_old * dt_factor
    dt_new = max(1e-6, min(dt_new, 1e-2))  # 限制在合理范围内
    
    return dt_new, dt_factor < 1.0

@njit
def calculate_fuel_heat_generation_numba(i, n_axial, dz):
    """计算燃料热生成 - Numba优化版本"""
    q_max = 1e6  # 最大线功率密度 (W/m)
    z = i * dz
    H = n_axial * dz
    # 使用余弦分布模拟轴向功率分布
    q = q_max * np.cos(np.pi * z / H)
    return q

@njit
def calculate_radial_heat_conduction_numba(T_c, T_f, dr):
    """计算径向热传导 - Numba优化版本"""
    k = 30.0  # 假设的导热系数 (W/m·K)
    # 简化的径向热传导计算
    q_rad = k * (T_c - T_f) / dr
    return q_rad

def set_boundary_conditions(self):
    """
    设置动态边界条件
    
    包括:
    1. 入口条件的动态调整
    2. 出口条件的处理
    3. 壁面条件的设定
    """
    try:
        # === 1. 入口边界条件 ===
        for j in range(self.n_channels):
            channel_type = self.channel_types[j+1]
            flow_area = self.channel_geometry[channel_type]['flow_area']
            
            # 1.1 入口空泡率 - 考虑启动过程
            if self.current_time < 0.1:  # 启动阶段
                self.void_fraction[0,j] = 0.0
            else:
                # 根据上游条件动态调整
                self.void_fraction[0,j] = min(0.1, self.void_fraction[1,j])
            
            # 1.2 入口速度分布 - 考虑子通道类型
            props = self.get_fluid_properties(self.inlet_pressure, self.inlet_temp, 'liquid')
            if props is not None:
                # 根据子通道类型调整质量流量
                if channel_type == 'corner':
                    factor = 0.7  # 角通道流量系数
                elif channel_type == 'edge':
                    factor = 0.9  # 边通道流量系数
                else:  # center
                    factor = 1.0  # 中心通道流量系数
                    
                # 计算入口速度
                self.liquid_velocity[0,j] = factor * self.inlet_mass_flux / props['density']
                self.vapor_velocity[0,j] = 1.2 * self.liquid_velocity[0,j]  # 气相速度略大
                
            # 1.3 入口温度和压力
            self.temperature[0,j] = self.inlet_temp
            self.pressure[0,j] = self.inlet_pressure
            
            # 1.4 入口焓值
            if props is not None:
                self.h_l[0,j] = props['enthalpy']
                props_v = self.get_fluid_properties(self.inlet_pressure, self.inlet_temp, 'vapor')
                if props_v is not None:
                    self.h_v[0,j] = props_v['enthalpy']
            
        # === 2. 出口边界条件 ===
        for j in range(self.n_channels):
            # 2.1 零梯度条件
            self.void_fraction[-1,j] = self.void_fraction[-2,j]
            self.liquid_velocity[-1,j] = self.liquid_velocity[-2,j]
            self.vapor_velocity[-1,j] = self.vapor_velocity[-2,j]
            
            # 2.2 出口压力指定
            self.pressure[-1,j] = self.outlet_pressure
            
            # 2.3 温度外推
            if self.n_axial > 2:
                self.temperature[-1,j] = 2*self.temperature[-2,j] - self.temperature[-3,j]
            else:
                self.temperature[-1,j] = self.temperature[-2,j]
            
            # 2.4 出口焓值外推
            if self.n_axial > 2:
                self.h_l[-1,j] = 2*self.h_l[-2,j] - self.h_l[-3,j]
                self.h_v[-1,j] = 2*self.h_v[-2,j] - self.h_v[-3,j]
            else:
                self.h_l[-1,j] = self.h_l[-2,j]
                self.h_v[-1,j] = self.h_v[-2,j]
            
        # === 3. 壁面边界条件 ===
        for i in range(self.n_axial):
            for j in range(self.n_channels):
                channel_type = self.channel_types[j+1]
                
                # 3.1 壁面无滑移条件
                if channel_type == 'corner':
                    self.liquid_velocity[i,j] *= 0.9  # 考虑强壁面效应
                    self.vapor_velocity[i,j] *= 0.9
                    
                # 3.2 壁面温度条件
                if hasattr(self, 'wall_temperature'):
                    self.T_c[i,j] = self.wall_temperature
                else:
                    # 使用简化的壁温计算
                    self.T_c[i,j] = self.temperature[i,j] + 30  # 假设温差30K
                    
                # 3.3 相变率边界条件
                if i == 0 or i == self.n_axial-1:
                    self.gamma[i,j] = 0.0
                    
        # === 4. 检查物理约束 ===
        self.void_fraction = np.clip(self.void_fraction, 0.0, 0.9999)
        self.liquid_velocity = np.maximum(self.liquid_velocity, 0.0)
        self.vapor_velocity = np.maximum(self.vapor_velocity, 0.0)
        self.pressure = np.maximum(self.pressure, 1e5)  # 最小压力1bar
        self.temperature = np.clip(self.temperature, 273.15, 647.0)  # 水的物理范围
        
    except Exception as e:
        logging.error(f"设置边界条件时出错: {str(e)}")
        logging.error(f"错误位置: {traceback.format_exc()}")

def apply_numerical_stabilization(self):
        """
        应用数值稳定化技术
        
        包括:
        1. 自适应时间步长控制
        2. 数值通量限制器
        3. 压力-速度耦合处理
        4. 数值振荡抑制
        """
        try:
            # 1. 自适应时间步长控制
            dt_cfl = self.calculate_cfl_timestep()
            dt_diff = self.calculate_diffusion_timestep()
            self.dt = min(dt_cfl, dt_diff)
            
            # 2. 应��通量��制器
            self.apply_flux_limiter()
            
            # 3. 压力-速度耦合处理
            self.apply_pressure_velocity_coupling()
            
            # 4. 数值振荡抑制
            self.apply_oscillation_damping()
            
        except Exception as e:
            logging.error(f"应用数值稳定化时出错: {str(e)}")
            
def calculate_cfl_timestep(self):
        """计算基于CFL条件的时间步长"""
        try:
            # 计算最大速度
            max_velocity = max(
                np.max(np.abs(self.liquid_velocity)),
                 np.max(np.abs(self.vapor_velocity))
            )
            
            
            # 考虑声速的影响
            props_l = self.get_fluid_properties(self.inlet_pressure, self.inlet_temp, 'liquid')
            if props_l is not None:
                sound_speed = np.sqrt(1.4 * self.inlet_pressure / props_l['density'])
                max_velocity = max(max_velocity, sound_speed)
                
            # 计算CFL时间步长
            dt_cfl = 0.5 * self.dz / max_velocity  # CFL数取0.5
            
            return dt_cfl
            
        except Exception as e:
            logging.error(f"计算CFL时间步长时出错: {str(e)}")
            return 1e-4  # 返回一个安全的默认值
            
def calculate_diffusion_timestep(self):
        """计算扩散稳定性时间步长"""
        try:
            # 获取流体物性
            props = self.get_fluid_properties(self.inlet_pressure, self.inlet_temp, 'liquid')
            if props is None:
                return 1e-4
                
            # 计算热扩散系数
            alpha = props.get('thermal_conductivity', 0.6) / \
                    (props['density'] * props.get('specific_heat', 4200))
                    
            # 计算动力粘性扩散系数
            nu = props['viscosity'] / props['density']
            
            # 取最扩散系数
            max_diff = max(alpha, nu)
            
            # 计算扩散时间步长
            dt_diff = 0.25 * self.dz**2 / max_diff  # 因子0.25确保稳定性
            
            return dt_diff
            
        except Exception as e:
            logging.error(f"计算扩散时间步长时出错: {str(e)}")
            return 1e-4
            
def apply_flux_limiter(self):
        """应用通量限制器以抑制数值振荡"""
        try:
            for i in range(1, self.n_axial-1):
                for j in range(self.n_channels):
                    # 计算空泡率梯度
                    r_void = (self.void_fraction[i,j] - self.void_fraction[i-1,j]) / \
                            (self.void_fraction[i+1,j] - self.void_fraction[i,j] + 1e-10)
                    
                    # 使用van Leer限制器
                    phi = (r_void + abs(r_void)) / (1 + abs(r_void))
                    
                    # 限制空泡率变化
                    self.void_fraction[i,j] = self.void_fraction[i-1,j] + \
                                            phi * (self.void_fraction[i+1,j] - self.void_fraction[i,j])
                                            
                    # 对速度场应用类似处理
                    r_vl = (self.liquid_velocity[i,j] - self.liquid_velocity[i-1,j]) / \
                           (self.liquid_velocity[i+1,j] - self.liquid_velocity[i,j] + 1e-10)
                    phi_vl = (r_vl + abs(r_vl)) / (1 + abs(r_vl))
                    self.liquid_velocity[i,j] = self.liquid_velocity[i-1,j] + \
                                              phi_vl * (self.liquid_velocity[i+1,j] - self.liquid_velocity[i,j])
                                              
        except Exception as e:
            logging.error(f"应用通量限制器时出错: {str(e)}")
            
def apply_pressure_velocity_coupling(self):
        """应用SIMPLE算法处理压力-速度耦合"""
        try:
            # 1. 求解动量预测方程
            u_star = np.copy(self.liquid_velocity)
            v_star = np.copy(self.vapor_velocity)
            
            # 2. 求解压力校正方程
            p_correction = self.solve_pressure_correction()
            
            # 3. 校正压力场
            self.pressure += 0.8 * p_correction  # 松弛因子0.8
            
            # 4. 校正速度场
            for i in range(1, self.n_axial-1):
                for j in range(self.n_channels):
                    # 速度校正
                    dp_dx = (p_correction[i+1,j] - p_correction[i-1,j]) / (2*self.dz)
                    props = self.get_fluid_properties(self.pressure[i,j], 
                                                    self.temperature[i,j], 'liquid')
                    if props is not None:
                        self.liquid_velocity[i,j] = u_star[i,j] - \
                                                  self.dt * dp_dx / props['density']
                        self.vapor_velocity[i,j] = v_star[i,j] - \
                                                 self.dt * dp_dx / (props['density'] * 0.1)
                                                 
        except Exception as e:
            logging.error(f"应用压力-速度耦合时出错: {str(e)}")
            
def apply_oscillation_damping(self):
        """应用数值振荡抑制技术"""
        try:
            # 1. 人工粘性项
            epsilon = 0.1  # 人工粘性系数
            
            for i in range(1, self.n_axial-1):
                for j in range(self.n_channels):
                    # 对空泡率应用人工粘性
                    d2alpha = (self.void_fraction[i+1,j] - 2*self.void_fraction[i,j] + 
                              self.void_fraction[i-1,j])
                    self.void_fraction[i,j] += epsilon * d2alpha
                    
                    # 对速度场应用人工粘性
                    d2vl = (self.liquid_velocity[i+1,j] - 2*self.liquid_velocity[i,j] + 
                            self.liquid_velocity[i-1,j])
                    self.liquid_velocity[i,j] += epsilon * d2vl
                    
                    d2vv = (self.vapor_velocity[i+1,j] - 2*self.vapor_velocity[i,j] + 
                            self.vapor_velocity[i-1,j])
                    self.vapor_velocity[i,j] += epsilon * d2vv
                    
            # 2. 限制物理范围
            self.void_fraction = np.clip(self.void_fraction, 0, 0.9999)
            self.liquid_velocity = np.maximum(self.liquid_velocity, 0)
            self.vapor_velocity = np.maximum(self.vapor_velocity, 0)
            self.pressure = np.maximum(self.pressure, 1e5)  # 最小压力1bar
            self.temperature = np.clip(self.temperature, 273.15, 647)  # 水的物理范围
            
        except Exception as e:
            logging.error(f"应用振荡抑制时出错: {str(e)}")

def solve_pressure_correction(self):
        """求解压力校正方程"""
        try:
            n = self.n_axial * self.n_channels
            A = lil_matrix((n, n))
            b = np.zeros(n)
            
            for i in range(1, self.n_axial-1):
                for j in range(self.n_channels):
                    idx = i * self.n_channels + j
                    
                    # 获取流体物性
                    props = self.get_fluid_properties(self.pressure[i,j], 
                                                    self.temperature[i,j], 'liquid')
                    if props is None:
                        continue
                        
                    # 构建系数矩阵
                    A[idx, idx] = -2
                    A[idx, idx-self.n_channels] = 1
                    A[idx, idx+self.n_channels] = 1
                    
                    # 构建源项（质量不平衡）
                    dm_in = props['density'] * self.liquid_velocity[i-1,j]
                    dm_out = props['density'] * self.liquid_velocity[i+1,j]
                    b[idx] = dm_in - dm_out
                    
            # 求解线性方程组
            p_correction = spsolve(A.tocsr(), b)
            
            # 重塑为二维数组
            return p_correction.reshape((self.n_axial, self.n_channels))
            
        except Exception as e:
            logging.error(f"求解压力校正方程时出错: {str(e)}")
            return np.zeros((self.n_axial, self.n_channels))

def calculate_heat_transfer(self, i, j):
        """
        计算综合传热过程
        
        包括:
        1. 轴向热传导
        2. 径向热传导
        3. 壁面传热
        4. 相间传热
        """
        try:
            # 获取局部条件
            P = self.pressure[i,j]
            T_l = self.T_l[i,j]
            T_v = self.T_v[i,j]
            T_w = self.T_c[i,j]
            
            # 1. 轴向热传导
            q_axial = self.calculate_axial_heat_conduction(i, j)
            
            # 2. 径向热传导
            q_radial = self.calculate_radial_heat_conduction(i, j)
            
            # 3. 壁面传热
            q_wall = self.calculate_wall_heat_transfer(i, j)
            
            # 4. 相间传热
            q_interface = self.calculate_interfacial_heat_transfer(i, j)
            
            return q_axial, q_radial, q_wall, q_interface
            
        except Exception as e:
            logging.error(f"计算综合传热时出错: {str(e)}")
            return 0.0, 0.0, 0.0, 0.0

def calculate_axial_heat_conduction(self, i, j):
        """计算轴向热传导"""
        try:
            if i == 0 or i == self.n_axial-1:
                return 0.0
                
            # 获取流体物性
            props_l = self.get_fluid_properties(self.pressure[i,j], self.T_l[i,j], 'liquid')
            if props_l is None:
                return 0.0
                
            # 计算温度梯度
            dT_dz = (self.T_l[i+1,j] - self.T_l[i-1,j]) / (2*self.dz)
            d2T_dz2 = (self.T_l[i+1,j] - 2*self.T_l[i,j] + self.T_l[i-1,j]) / (self.dz**2)
            
            # 计算轴向热流密度
            q_axial = -props_l['thermal_conductivity'] * d2T_dz2
            
            return q_axial
            
        except Exception as e:
            logging.error(f"计算轴向热传导时出错: {str(e)}")
            return 0.0

def calculate_radial_heat_conduction(self, i, j):
        """计算径向热传导"""
        try:
            # 获取子通道类型和几何参数
            channel_type = self.channel_types[j+1]
            flow_area = self.channel_geometry[channel_type]['flow_area']
            wetted_perimeter = self.channel_geometry[channel_type]['wetted_perimeter']
            
            # 计算特征长度
            L_char = 4 * flow_area / wetted_perimeter
            
            # 获取流体物性
            props_l = self.get_fluid_properties(self.pressure[i,j], self.T_l[i,j], 'liquid')
            if props_l is None:
                return 0.0
                
            # 计算径向温度梯度
            dT_dr = (self.T_c[i,j] - self.T_l[i,j]) / (L_char/2)
            
            # 计算径向热流密度
            q_radial = -props_l['thermal_conductivity'] * dT_dr
            
            return q_radial
            
        except Exception as e:
            logging.error(f"计算径向热传导时出错: {str(e)}")
            return 0.0

def calculate_wall_heat_transfer(self, i, j):
        """计算壁面传热"""
        try:
            # 获取局部条件
            T_w = self.T_c[i,j]
            T_l = self.T_l[i,j]
            T_v = self.T_v[i,j]
            P = self.pressure[i,j]
            
            # 计算饱和温度
            T_sat = PropsSI('T', 'P', P, 'Q', 0, 'Water')
            
            # 根据传热机制选择合适的换热系数
            if T_l < T_sat and T_w > T_sat:
                # 过冷沸腾
                h = self.calculate_subcooled_boiling_htc(T_w, T_sat, T_l, P)
            elif T_l >= T_sat and T_w > T_sat:
                # 饱和核态沸腾
                h = self.calculate_nucleate_boiling_htc(T_w, T_sat, P)
            else:
                # 单相对流换热
                h = self.calculate_single_phase_htc(T_w, P, 'liquid')
                
            # 计算壁面热流密度
            q_wall = h * (T_w - T_l)
            
            return q_wall
            
        except Exception as e:
            logging.error(f"计算壁面传热时出错: {str(e)}")
            return 0.0

def calculate_interfacial_heat_transfer(self, i, j):
        """计算相间传热"""
        try:
            # 获取局部条件
            T_l = self.T_l[i,j]
            T_v = self.T_v[i,j]
            P = self.pressure[i,j]
            alpha = self.void_fraction[i,j]
            
            # 获取流体物性
            props_l = self.get_fluid_properties(P, T_l, 'liquid')
            props_v = self.get_fluid_properties(P, T_v, 'vapor')
            
            if props_l is None or props_v is None:
                return 0.0
                
            # 计算界面面积密度
            d_bubble = 0.001  # 假设气泡直径为1mm
            a_i = 6 * alpha / d_bubble
            
            # 计算相对速度
            v_r = abs(self.vapor_velocity[i,j] - self.liquid_velocity[i,j])
            
            # 计算界面换热系数
            Re_i = props_l['density'] * v_r * d_bubble / props_l['viscosity']
            Pr_l = props_l['viscosity'] * props_l['specific_heat'] / props_l['thermal_conductivity']
            Nu_i = 2 + 0.6 * Re_i**0.5 * Pr_l**0.33
            h_i = Nu_i * props_l['thermal_conductivity'] / d_bubble
            
            # 计算相间传热量
            q_interface = h_i * a_i * (T_v - T_l)
            
            return q_interface
            
        except Exception as e:
            logging.error(f"计算相间传热时出错: {str(e)}")
            return 0.0

def update_temperature_fields(self):
        """更新温度场"""
        try:
            for i in range(1, self.n_axial-1):
                for j in range(self.n_channels):
                    # 获取各传热分量
                    q_axial, q_radial, q_wall, q_interface = self.calculate_heat_transfer(i, j)
                    
                    # 获取流体物性
                    props_l = self.get_fluid_properties(self.pressure[i,j], self.T_l[i,j], 'liquid')
                    props_v = self.get_fluid_properties(self.pressure[i,j], self.T_v[i,j], 'vapor')
                    
                    if props_l is None or props_v is None:
                        continue
                        
                    # 计算温度变化率
                    dTl_dt = (q_axial + q_radial + q_wall + q_interface) / \
                             (props_l['density'] * props_l['specific_heat'])
                         
                    dTv_dt = -q_interface / (props_v['density'] * props_v['specific_heat'])
                    
                    # 更新温度
                    self.T_l[i,j] += self.dt * dTl_dt
                    self.T_v[i,j] += self.dt * dTv_dt
                    
                    # 更新壁温
                    self.update_wall_temperature(i, j)
                    
        except Exception as e:
            logging.error(f"更新温度场时出错: {str(e)}")

def update_wall_temperature(self, i, j):
        """更新壁面温度"""
        try:
            # 获取子通道类型和几何参数
            channel_type = self.channel_types[j+1]
            heated_perimeter = self.channel_geometry[channel_type]['heated_perimeter']
            
            # 计算燃料产热
            q_fuel = self.calculate_fuel_heat_generation(i, j)
            
            # 计算壁面导热
            k_wall = 16.0  # 假设壁面材料为不锈钢
            delta_wall = 0.001  # 假设壁厚1mm
            
            # 计算新的壁温
            self.T_c[i,j] = self.T_l[i,j] + \
                            (q_fuel * heated_perimeter * delta_wall) / k_wall
                        
        except Exception as e:
            logging.error(f"更新壁面温度时出错: {str(e)}")

def calculate_fuel_heat_generation(self, i, j):
        """计算燃料产热"""
        try:
            # 使用余弦分布模拟轴向功率分布
            z = i * self.dz
            q_max = 1e6  # 最大线功率密度
            
            # 考虑子通道类型的功率因子
            channel_type = self.channel_types[j+1]
            if channel_type == 'corner':
                power_factor = 0.8
            elif channel_type == 'edge':
                power_factor = 0.9
            else:  # center
                power_factor = 1.0
                
            # 计算局部功率密度
            q = q_max * power_factor * np.cos(np.pi * z / self.channel_length)
            
            return q
            
        except Exception as e:
            logging.error(f"计算燃料产热时出错: {str(e)}")
            return 0.0

@njit(parallel=True)
def compute_grid_points_numba(n_axial, n_channels,
                            void_fraction, prev_void_fraction,
                            vapor_velocity, liquid_velocity,
                            prev_liquid_velocity, prev_vapor_velocity,
                            pressure, temperature,
                            time_step, dz):
    """使用Numba优化的网格点计算"""
    new_void = np.zeros_like(void_fraction)
    new_vl = np.zeros_like(liquid_velocity)
    new_vv = np.zeros_like(vapor_velocity)
    
    for i in prange(1, n_axial-1):
        for j in prange(n_channels):
            # 计算空间导数
            d_alpha_dx = (void_fraction[i+1,j] - void_fraction[i-1,j]) / (2*dz)
            d_vl_dx = (liquid_velocity[i+1,j] - liquid_velocity[i-1,j]) / (2*dz)
            d_vv_dx = (vapor_velocity[i+1,j] - vapor_velocity[i-1,j]) / (2*dz)
            
            # 计算压力梯度
            dp_dx = (pressure[i+1,j] - pressure[i-1,j]) / (2*dz)
            
            # 更新空泡率
            new_void[i,j] = prev_void_fraction[i,j] - time_step * (
                vapor_velocity[i,j] * d_alpha_dx
            )
            
            # 更新液相速度
            new_vl[i,j] = prev_liquid_velocity[i,j] - time_step * (
                liquid_velocity[i,j] * d_vl_dx +
                dp_dx / 1000.0 +  # 假设液相密度约为1000 kg/m³
                9.81  # 重力项
            )
            
            # 更新气相速度
            new_vv[i,j] = prev_vapor_velocity[i,j] - time_step * (
                vapor_velocity[i,j] * d_vv_dx +
                dp_dx / 1.0 +  # 假设气相密度约为1 kg/m³
                9.81  # 重力项
            )
            
            # 限制物理范围
            new_void[i,j] = min(max(new_void[i,j], 0.0), 0.9999)
            new_vl[i,j] = max(new_vl[i,j], 0.0)
            new_vv[i,j] = max(new_vv[i,j], 0.0)
            
    # 处理边界条件
    new_void[0,:] = void_fraction[0,:]
    new_void[-1,:] = void_fraction[-2,:]
    new_vl[0,:] = liquid_velocity[0,:]
    new_vl[-1,:] = liquid_velocity[-2,:]
    new_vv[0,:] = vapor_velocity[0,:]
    new_vv[-1,:] = vapor_velocity[-2,:]
    
    return new_void, new_vl, new_vv

if __name__ == "__main__":
    # 配置日志
    logging.basicConfig(
        level=logging.WARNING,
        format='%(asctime)s - %(name)s - %(levelname)s - %(message)s',
        force=True  # 强制重新配置
    )

    # 创建文件处理器
    file_handler = logging.FileHandler('subchannel_analysis.log', mode='w', encoding='utf-8')
    file_handler.setLevel(logging.WARNING)
    file_handler.setFormatter(logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s'))

    # 创建控制台处理器
    console_handler = logging.StreamHandler()
    console_handler.setLevel(logging.WARNING)
    console_handler.setFormatter(logging.Formatter('%(levelname)s: %(message)s'))

    # 获取根日志记录器
    logger = logging.getLogger()
    logger.handlers = []  # 清除所有处理器
    logger.addHandler(file_handler)
    logger.addHandler(console_handler)

    # 设置CoolProp警告过滤
    logging.getLogger('CoolProp').setLevel(logging.ERROR)
    
    try:
        # 创建分析实例
        analysis = SubchannelAnalysis(n_axial=20)
        
        # 运行瞬态模拟
        analysis.run_transient_simulation()
        
        # 显示瞬态结果
        analysis.plot_transient_results()
        
        # 显示最终状态的详细分布
        analysis.plot_results()
        
    except Exception as e:
        logging.error(f"程序执行出错: {str(e)}")
        logging.error(f"错误位置: {traceback.format_exc()}")